Polysubstance and tr. completion

Analyze the association between polysubstance at admission and tr. compleiton longitudinally along treatments, accounting for irregular observations.

Author

Andrés González Santa Cruz

Published

February 16, 2025

Setup

Code
rm(list = ls()) 
folder_path <- ifelse(dir.exists("E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/"),
                      "E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/",
                      "G:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/")
load(paste0(folder_path,"an_grant_23_24_3.RData"))

Packages

Code
system("fc-cache -f -v")
[1] 0
Code
if(!require(pacman)){install.packages("pacman");library(pacman)}
pacman::p_unlock(lib.loc = .libPaths()) #para no tener problemas reinstalando paquetes


if (getRversion() != "4.1.2") { stop("Requires R version 4.1.2. Actual: ", getRversion()) }


if(!require(geepack)){install.packages("geepack");library(geepack)}
if(!require(geeM)){install.packages("geeM");library(geeM)}
if(!require(tidyverse)){install.packages("tidyverse");library(tidyverse)}
if(!require(MASS)){install.packages("MASS");library(MASS)}
if(!require(geeasy)){install.packages("geeasy");library(geeasy)}
if(!require(MuMIn)){remotes::install_version("MuMIn", "1.46.0");library(MuMIn)}
if(!require(tableone)){install.packages("tableone");library(tableone)}
if(!require(knitr)){install.packages("knitr");library(knitr)}
if(!require(emmeans)){install.packages("emmeans");library(emmeans)}
if(!require(biostat3)){install.packages("biostat3");library(biostat3)}
if(!require(rms)){install.packages("rms");library(rms)}
if(!require(meta)){install.packages("meta");library(meta)}
try(if(!require(dmetar)){install.packages("dmetar");library(dmetar)})
Error : package 'dmetar' is not available
Code
if(!require(metafor)){install.packages("metafor");library(metafor)}
try(if(!require(extrafont)){install.packages("extrafont");library(extrafont)})
try(if(!require(showtext)){install.packages("showtext");library(showtext)})
try(if(!require(mice)){install.packages("mice");library(mice)})
try(if(!require(htmltools)){install.packages("htmltools");library(htmltools)})
if(!require(chisq.posthoc.test)){devtools::install_github("ebbertd/chisq.posthoc.test")}
if(!require(DescTools)){install.packages("DescTools")}
if(!require(vcd)){install.packages("vcd")}

try(extrafont::font_import(pattern = "Times New Roman", prompt = FALSE))
Error in data.frame(fontfile = ttfiles, FontName = "", stringsAsFactors = FALSE) : 
  arguments imply differing number of rows: 0, 1
Code
extrafont::loadfonts(device="win")
# Agregar la ruta a la fuente "Times New Roman"
font_add("Times New Roman", "C:/Windows/Fonts/times.ttf")

# Habilitar showtext para gráficos
showtext::showtext_auto()

#https://gist.github.com/avallecam/56af06f46e5544c3af0f46344df20989 
#https://stackoverflow.com/questions/49089476/any-updates-to-model-negative-binomial-distribution-data-with-gee-in-r
#https://www.rdocumentation.org/packages/geepack/versions/1.3.4/topics/geese
#https://stackoverflow.com/questions/13946540/negative-binomial-in-gee?rq=1


#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

summary2.geem <- function(object, exponentiate = FALSE, digits = 2, ...) {
  if (!is.list(object) || !"beta" %in% names(object)) {
    stop("Invalid object: Expected a model object with beta coefficients.")
  }
  
  # Initialize coefficients matrix
  Coefs <- matrix(NA, nrow = length(object$beta), ncol = 8)
  rownames(Coefs) <- object$coefnames
  colnames(Coefs) <- c("Term", "Estimates", "Model SE", "Robust SE", "Wald", "CI.lower", "CI.upper", "p.value")
  
  Coefs[, "Estimates"] <- object$beta
  
  # Handling the variance calculations
  naive <- is.null(object$var) || is.character(object$var) || any(diag(object$var) < 0)
  if (naive) {
    warning("Robust variance estimate not available or invalid. Using model-based SE.")
    Coefs[, "Robust SE"] <- sqrt(diag(object$naiv.var))
  } else {
    Coefs[, "Robust SE"] <- sqrt(diag(object$var))
  }
  
  Coefs[, "Model SE"] <- sqrt(diag(object$naiv.var))  # Model-based SE
  
  # Calculate Wald statistics
  Coefs[, "Wald"] <- Coefs[, "Estimates"] / Coefs[, "Robust SE"]
  
  # 95% CI calculations
  z_value <- qnorm((1 + 0.95) / 2)
  Coefs[, "CI.lower"] <- Coefs[, "Estimates"] - z_value * Coefs[, "Robust SE"]
  Coefs[, "CI.upper"] <- Coefs[, "Estimates"] + z_value * Coefs[, "Robust SE"]
  
  # P-value calculation
  Coefs[, "p.value"] <- round(2 * pnorm(-abs(Coefs[, "Wald"])), digits = 4)
  
  # Exponentiate estimates and CIs if specified
  if (exponentiate) {
    Coefs[, c("Estimates", "CI.lower", "CI.upper")] <- exp(Coefs[, c("Estimates", "CI.lower", "CI.upper")])
  }
  
  # Round numerical output
  decimal_format <- sprintf("%%1.%df", digits)
  
  Coefs[, c("Estimates", "Model SE", "Robust SE", "Wald", "CI.lower", "CI.upper")] <- round(Coefs[, c("Estimates", "Model SE", "Robust SE", "Wald", "CI.lower", "CI.upper")], digits = digits)
  Coefs[, c("Estimates", "Model SE", "Robust SE", "Wald", "CI.lower", "CI.upper")] <- sapply(Coefs[, c("Estimates", "Model SE", "Robust SE", "Wald", "CI.lower", "CI.upper")], function(col) {
    as.character(sprintf(decimal_format, as.numeric(col)))
  })
  # Convert matrix to data frame and keep row names
  result_df <- as.data.frame(Coefs)
  result_df[, "Term"] <- rownames(Coefs)
  row.names(result_df) <- NULL
  return(result_df)
}

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

tidy_geese_model <- function(geese_model, conf.level = 0.95) {
  # Extract coefficients and their names
  coefficients <- geese_model$beta
  coefficient_names <- names(coefficients)
  
  # Extract variance of beta coefficients and compute standard errors
  variance_beta <- geese_model$vbeta
  std_errors <- sqrt(diag(variance_beta))
  
  # Calculate z-values and p-values
  z_values <- coefficients / std_errors
  p_values <- 2 * pnorm(-abs(z_values), lower.tail = TRUE)
  
  # Calculate confidence intervals
  alpha <- 1 - conf.level
  z <- qnorm(1 - alpha / 2)
  lower_ci <- coefficients - z * std_errors
  upper_ci <- coefficients + z * std_errors
  
  # Exponentiate coefficients and confidence intervals for odds ratios
  exp_coefficients <- exp(coefficients)
  exp_lower_ci <- exp(lower_ci)
  exp_upper_ci <- exp(upper_ci)
  
  # Create a tidy data frame
  tidy_data <- tibble::tibble(
    Term = coefficient_names,
    Estimate = exp_coefficients,
    Std.Error = std_errors,
    z.value = z_values,
    p.value = p_values,
    CI.Lower = exp_lower_ci,
    CI.Upper = exp_upper_ci
  ) %>% dplyr::mutate_at(c("Estimate", "Std.Error", "z.value", "CI.Lower","CI.Upper"),~sprintf("%1.2f",.))%>% dplyr::mutate_at(c("p.value"),~sprintf("%1.4f",.))
  
  return(tidy_data)
}

options(knitr.kable.NA = '')


#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

knitr::knit_hooks$set(time_it = local({
  now <- NULL
  function(before, options) {
    if (before) {
      # record the current time before each chunk
      now <<- Sys.time()
    } else {
      # calculate the time difference after a chunk
      res <- ifelse(difftime(Sys.time(), now)>(60^2),difftime(Sys.time(), now)/(60^2),difftime(Sys.time(), now)/(60^1))
      # return a character string to show the time
      x<-ifelse(difftime(Sys.time(), now)>(60^2),paste("Time for this code chunk to run:", round(res,1), "hours"),paste("Time for this code chunk to run:", round(res,1), "minutes"))
      paste('<div class="message">', gsub('##', '\n', x),'</div>', sep = '\n')
    }
  }
}))
knitr::opts_chunk$set(time_it = TRUE)

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

format_cells <- function(df, rows ,cols, value = c("italics", "bold", "strikethrough")){
  
  # select the correct markup
  # one * for italics, two ** for bold
  map <- setNames(c("*", "**", "~~"), c("italics", "bold", "strikethrough"))
  markup <- map[value]  
  
  for (r in rows){
    for(c in cols){
      
      # Make sure values are not factors
      df[[c]] <- as.character( df[[c]])
      
      # Update formatting
      df[r, c] <- ifelse(nchar(df[r, c])==0,"",paste0(markup, gsub(" ", "", df[r, c]), markup))
    }
  }
  
  return(df)
}
#To produce line breaks in messages and warnings
knitr::knit_hooks$set(
  error = function(x, options) {
    paste('\n\n<div class="alert alert-danger">',
          gsub('##', '\n', gsub('^##\ Error', '**Error**', x)),
          '</div>', sep = '\n')
  },
  warning = function(x, options) {
    paste('\n\n<div class="alert alert-warning">',
          gsub('##', '\n', gsub('^##\ Warning:', '**Warning**', x)),
          '</div>', sep = '\n')
  },
  message = function(x, options) {
    paste('<div class="message">',
          gsub('##', '\n', x),
          '</div>', sep = '\n')
  }
)

#_#_#_#_#_#_#_#_#_#_#_#_#_
invisible("Function to format CreateTableOne into a database")

as.data.frame.TableOne <- function(x, ...) {capture.output(print(x,
                                                                 showAllLevels = TRUE, varLabels = T,...) -> x)
  y <- as.data.frame(x)
  y$characteristic <- dplyr::na_if(rownames(x), "")
  y <- y %>%
    fill(characteristic, .direction = "down") %>%
    dplyr::select(characteristic, everything())
  rownames(y) <- NULL
  y}
#_#_#_#_#_#_#_#_#_#_#_#_#_
# Austin, P. C. (2009). The Relative Ability of Different Propensity 
# Score Methods to Balance Measured Covariates Between 
# Treated and Untreated Subjects in Observational Studies. Medical 
# Decision Making. https://doi.org/10.1177/0272989X09341755
smd_bin <- function(x,y){
  z <- x*(1-x)
  t <- y*(1-y)
  k <- sum(z,t)
  l <- k/2
  
  return((x-y)/sqrt(l))
  
}

#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_

counting_process_output<-
  function(object=NULL){
    broom::tidy(object$m, exponentiate=T, conf.int=T) %>% 
      dplyr::mutate(across(c("estimate","std.error","robust.se","statistic","conf.low","conf.high"),~sprintf("%1.2f",.))) %>% 
      dplyr::mutate(p.value= sprintf("%1.4f",p.value))
  }

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:

tidy_coxph <- function(model, conf.level = 0.95) {
  # Extract the coefficients (a named numeric vector)
  coefs <- model$coefficients
  
  # Compute robust standard errors from the robust variance matrix
  robust_se <- sqrt(diag(model$var))
  
  # Compute z statistics and two-sided p-values
  z <- coefs / robust_se
  p_value <- 2 * (1 - pnorm(abs(z)))
  
  # Compute the critical z value for the specified confidence level
  z_crit <- qnorm(1 - (1 - conf.level) / 2)
  
  # Exponentiate coefficients to get hazard ratios
  hazard_ratio <- exp(coefs)
  
  # Compute the lower and upper confidence intervals (exponentiated)
  lower_ci <- exp(coefs - z_crit * robust_se)
  upper_ci <- exp(coefs + z_crit * robust_se)
  
  # Build and return a tidy data frame
  result <- data.frame(
    term         = names(coefs),
    coef         = coefs,
    robust_se    = robust_se,
    hazard_ratio = hazard_ratio,
    lower_ci     = lower_ci,
    upper_ci     = upper_ci,
    z            = z,
    p_value      = p_value,
    stringsAsFactors = FALSE
  )
  
  return(result)
}

format_hr_ci <- function(tidy_table, digits = 2) {
  # Build a format string dynamically using the desired number of digits.
  fmt <- paste0("%.", digits, "f (%.", digits, "f, %.", digits, "f)")
  
  # Create the new column using sprintf to format each row.
  tidy_table$hr_ci <- with(tidy_table, sprintf(fmt, hazard_ratio, lower_ci, upper_ci))
  
  return(tidy_table)
}

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:


if(.Platform$OS.type == "windows") withAutoprint({
  memory.size()
  memory.size(TRUE)
  memory.limit()
})
> memory.size()
[1] 7546.27
> memory.size(TRUE)
[1] 7546.88
> memory.limit()
[1] 32563
Code
memory.limit(size=56000)
[1] 56000
Code
invisible("At the beginning is in wide format")
invisible("At the beginning is in wide format")

nrow(Base_fiscalia_v15f_grant_23_24_long)
#[1] 109756
length(unique(Base_fiscalia_v15f_grant_23_24_long$hash_key))
#[1] 85048


Base_fiscalia_v15f_grant_23_24_long |> 
    dplyr::filter(motivodeegreso_mod_imp %in% c("En curso", "Muerte", "Derivación")) |>
    (\(df) {
        print(message(paste0("Dead, Referrals, ongoing treatments, Entries: ", nrow(df))))
        print(message(paste0("Dead, Referrals, ongoing treatments, HASHs: ", nrow(distinct(df,hash_key)))))
        
    })()

Base_fiscalia_v15f_grant_23_24_long |> 
    dplyr::filter(!motivodeegreso_mod_imp %in% c("En curso", "Muerte", "Derivación")) |>
    (\(df) {
        print(message(paste0("Dead, Referrals, ongoing treatments, Entries: ", nrow(df))))
        print(message(paste0("Dead, Referrals, ongoing treatments, HASHs: ", nrow(distinct(df,hash_key)))))
        
    })()

length(unlist(unique(Base_fiscalia_v15f_grant_23_24_long2$hash_key)))
 
#[exclude ongoing treatments and external referrals]
length(unique(Base_fiscalia_v15f_grant_23_24$hash_key))

length(unique(Base_fiscalia_v15f_grant_23_24_long2$hash_key))
#72404
nrow(Base_fiscalia_v15f_grant_23_24_long2)
#[1] 90075
#
Base_fiscalia_v15f_grant_23_24_long2 |> 
    dplyr::filter(!hash_key %in% Base_fiscalia_v15f_grant_23_24_long2_miss_proc_multtr$hash_key
) |>
    (\(df) {
        print(message(paste0("Only one treatment: ", nrow(df))))
        print(message(paste0("Only one treatment: ", nrow(distinct(df,hash_key)))))
        
    })()

Time for this code chunk to run: 0 minutes

Code
library(DiagrammeR)
gr<-
grViz("
digraph flowchart {
  graph [layout = dot, rankdir = TB]

  # General node styling
  node [fontname = Times, shape = rectangle, fontsize = 10, style = filled, fillcolor = transparent]

  # Main flow nodes
  original [label = 'Original C1 Dataset\\n(n = 163,146;\\nPatients = 85,722)', fillcolor = lightgray]
  c1_dataset [label = 'C1 Dataset\\n(n = 109,756;\\nPatients = 85,048)']
  after_discard [label = 'C1 Dataset\\n(n = 90,075;\\nPatients = 72,404)']
  final_dataset [label = 'Final C1 Dataset\\n(n = 30,988;\\nPatients = 13,317)', fillcolor = lightgray]

  # Discard nodes (aligned between main flow steps)
  discard_referrals [label = '&#8226;Discard external referrals or ongoing treatments\\l(n = 19,681; Patients = 18,450)\\l']
  discard_duplicates [label = '&#8226;Discard duplicated entries\\l&#8226;Intermediate events of treatment (continuous referrals)\\l']
  discard_single [label = '&#8226;Discard admitted to only one treatment(Patients = 59,087)\\l']

  # Invisible vertices for middle line
  v1 [shape = point, width = 0, style = invis]
  v2 [shape = point, width = 0, style = invis]
  v3 [shape = point, width = 0, style = invis]

  # Main flow edges (vertical line)
  original -> v1 [arrowhead = none]
  v1 -> c1_dataset
  c1_dataset -> v2 [arrowhead = none]
  v2 -> after_discard
  after_discard -> v3 [arrowhead = none]
  v3 -> final_dataset

  # Discard connections (from the middle line)
  v1 -> discard_referrals
  v2 -> discard_duplicates
  v3 -> discard_single

  # Alignment
  { rank = same; discard_referrals; v1 }
  { rank = same; discard_duplicates; v2 }
  { rank = same; discard_single; v3 }
}
", 
width = 800,
height = 600)
gr
Code
unlink(paste0(folder_path, "_doc/_flowchart_files"), recursive = TRUE)
htmlwidgets::saveWidget(gr, paste0(folder_path, "_doc/_flowchart.html"))
webshot::webshot(paste0(folder_path, "_doc/_flowchart.html"), 
                 paste0(folder_path, "_doc/_flowchart.png"),
                 vwidth = 300*1.2, vheight = 300,  zoom=10, expand=100)  # Prueba con diferentes coordenadas top, left, width, and height.

Time for this code chunk to run: 0.2 minutes

0.a. Descriptives

Code
invisible("2024-05-21: test different")
data_mine_miss_restr_proc2exp<-
data_mine_miss_restr_proc2 %>% 
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>% 
  {if (nrow(.) > nrow(data_mine_miss_restr_proc2)) {message("left join added rows"); .} else .} %>% 
  dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>% 
  dplyr::mutate(
    susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
    susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
    susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
    susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
    susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
    susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
    susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
    susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
  )

#"susinidumrec_otr", "susinidum_coc", "susinidum_pbc", "susinidum_mar",

  #janitor::tabyl(n_off_acq)

#(Post-Treatment)
# n_post_off_acq
# n_post_off_vio
# n_post_off_sud
# n_post_off_ot
# n_post_off

invisible("To account for variability by treatment setting, we stratified the analysis by setting: basic ambulatory GP intensive ambulatory GP residential WO intensive ambulatory WO residential")
table(subset(Base_fiscalia_v15f_grant_23_24_long2_miss_proc_multtr, is_first_occurrence==1, select= "tipo_de_plan_2_mod"))

       basic ambulatory GP intensive ambulatory          GP residential 
                   4360                    4998                    2178 
WO intensive ambulatory          WO residential 
                    745                    1036 
Code
data_mine_miss_restr_proc2exp %>% 
  group_by(hash_key) %>% 
  summarise(min(fech_ing_num)) %>% nrow()
[1] 13317
Code
#13317
subset(data_mine_miss_restr_proc2 , is_first_occurrence==1) %>% nrow()
[1] 13317
Code
#13317
data_mine_miss_restr_proc2_baseline<-
  subset(data_mine_miss_restr_proc2exp , is_first_occurrence==1) 


variables_vector <- c("tr_outcome", "policonsumo2", "comp_bpsc_y3_severe",
                      "less_90d_tr1", "log_dias_treat_imp_sin_na",
                      "edad_al_ing_1", "ano_nac_corr", "susinidumrec_coc", 
                      "susinidumrec_pbc", "susinidumrec_mar", "susinidumrec_otr", 
                      "psycom_dum_study", "psycom_dum_with", "freq_cons_dum_5day", 
                      "cond_oc_dum_2inact", "cond_oc_dum_3unemp", "susprindumrec_coc", 
                      "susprindumrec_pbc", "susprindumrec_mar", "susprindumrec_otr")#,"lag_comp_bpsc_y3_severe", "lag_less_90d_tr1", "log_lag_dias_treat_imp_sin_na", "lag_policonsumo2", "lag_tr_outcome")

Time for this code chunk to run: 0 minutes

Code
attr(data_mine_miss_restr_proc2_baseline$tr_outcome,"label") <- "Complete status of treatment (binary)"
attr(data_mine_miss_restr_proc2_baseline$tr_outcome,"label") <- "Complete status of treatment (binary)"
attr(data_mine_miss_restr_proc2_baseline$policonsumo2,"label") <- "Polysubstance use"
attr(data_mine_miss_restr_proc2_baseline$comp_bpsc_y3_severe,"label") <- "Biopsychosocial compromise (Severe)"
attr(data_mine_miss_restr_proc2_baseline$less_90d_tr1,"label") <- "Treatment duration (binary) (<90 days)"
attr(data_mine_miss_restr_proc2_baseline$log_dias_treat_imp_sin_na,"label") <- "Treatment duration (log-scaled days)"
attr(data_mine_miss_restr_proc2_baseline$edad_al_ing_1,"label") <- "Age at admission to treatment"
attr(data_mine_miss_restr_proc2_baseline$ano_nac_corr,"label") <- "Birth year"
attr(data_mine_miss_restr_proc2_baseline$susinidumrec_coc,"label") <- "Primary substance (initial diagnosis): cocaine hydrochloride"
attr(data_mine_miss_restr_proc2_baseline$susinidumrec_pbc,"label") <- "Primary substance (initial diagnosis): cocaine base paste"
attr(data_mine_miss_restr_proc2_baseline$susinidumrec_mar,"label") <- "Primary substance (initial diagnosis): marijuana"
attr(data_mine_miss_restr_proc2_baseline$susinidumrec_otr,"label") <- "Primary substance (initial diagnosis): other"
attr(data_mine_miss_restr_proc2_baseline$psycom_dum_study,"label") <- "Psychiatric comorbidity (ICD-10): In study"
attr(data_mine_miss_restr_proc2_baseline$psycom_dum_with,"label") <- "Psychiatric comorbidity (ICD-10): Diagnosed"
attr(data_mine_miss_restr_proc2_baseline$freq_cons_dum_5day,"label") <- "Daily frequence of primary substance use at admission"
attr(data_mine_miss_restr_proc2_baseline$cond_oc_dum_2inact,"label") <- "Occupational Status: Inactive"
attr(data_mine_miss_restr_proc2_baseline$cond_oc_dum_3unemp,"label") <- "Occupational Status: Unemployed"
attr(data_mine_miss_restr_proc2_baseline$susprindumrec_coc,"label") <- "Primary substance at admission to treatment: cocaine hydrochloride"
attr(data_mine_miss_restr_proc2_baseline$susprindumrec_pbc,"label") <- "Primary substance at admission to treatment: cocaine base paste"
attr(data_mine_miss_restr_proc2_baseline$susprindumrec_mar,"label") <- "Primary substance at admission to treatment: marijuana"
attr(data_mine_miss_restr_proc2_baseline$susprindumrec_otr,"label") <- "Primary substance at admission to treatment: other"
attr(data_mine_miss_restr_proc2_baseline$tipo_de_plan_2_mod,"label") <- "Treatment setting"

invisible("Table for imputed dataset with selected covariates")

library(tableone)

vars <- setdiff(c(variables_vector,"tipo_de_plan_2_mod"), "policonsumo2")
factorVars <- setdiff(variables_vector, c("policonsumo2", "edad_al_ing_1", "log_dias_treat_imp_sin_na", "ano_nac_corr"))

tableOne <- CreateTableOne(vars = vars, data = data_mine_miss_restr_proc2_baseline,
                           factorVars = factorVars, strata = "policonsumo2", 
                           addOverall = TRUE,  # <-- This is the key addition
                           includeNA = TRUE,
                           test = TRUE) # Remove statistical tests for speed

print(tableOne, smd = TRUE) # Calculate SMD separately if needed
                                       Stratified by policonsumo2
                                        Overall         0              
  n                                       13317            2383        
  tr_outcome = 1 (%)                      10448 (78.5)     1833 (76.9) 
  comp_bpsc_y3_severe = 1 (%)              5496 (41.3)      690 (29.0) 
  less_90d_tr1 = 1 (%)                     3269 (24.5)      567 (23.8) 
  log_dias_treat_imp_sin_na (mean (SD))    5.00 (1.07)     5.03 (0.93) 
  edad_al_ing_1 (mean (SD))               33.80 (9.46)    38.29 (11.32)
  ano_nac_corr (mean (SD))              1979.02 (9.54)  1975.34 (11.02)
  susinidumrec_coc = 1 (%)                  627 ( 4.7)      108 ( 4.5) 
  susinidumrec_pbc = 1 (%)                 1123 ( 8.4)      251 (10.5) 
  susinidumrec_mar = 1 (%)                 4348 (32.6)      483 (20.3) 
  susinidumrec_otr = 1 (%)                  295 ( 2.2)       57 ( 2.4) 
  psycom_dum_study = 1 (%)                 2653 (19.9)      420 (17.6) 
  psycom_dum_with = 1 (%)                  5836 (43.8)      986 (41.4) 
  freq_cons_dum_5day = 1 (%)               6242 (46.9)     1013 (42.5) 
  cond_oc_dum_2inact = 1 (%)               2410 (18.1)      468 (19.6) 
  cond_oc_dum_3unemp = 1 (%)               5271 (39.6)      764 (32.1) 
  susprindumrec_coc = 1 (%)                2370 (17.8)      292 (12.3) 
  susprindumrec_pbc = 1 (%)                6898 (51.8)      902 (37.9) 
  susprindumrec_mar = 1 (%)                 717 ( 5.4)       64 ( 2.7) 
  susprindumrec_otr = 1 (%)                 188 ( 1.4)       43 ( 1.8) 
  tipo_de_plan_2_mod (%)                                               
     basic ambulatory                      4360 (32.7)     1040 (43.6) 
     GP intensive ambulatory               4998 (37.5)      786 (33.0) 
     GP residential                        2178 (16.4)      272 (11.4) 
     WO intensive ambulatory                745 ( 5.6)      138 ( 5.8) 
     WO residential                        1036 ( 7.8)      147 ( 6.2) 
                                       Stratified by policonsumo2
                                        1               p      test SMD   
  n                                       10934                           
  tr_outcome = 1 (%)                       8615 (78.8)   0.047       0.045
  comp_bpsc_y3_severe = 1 (%)              4806 (44.0)  <0.001       0.315
  less_90d_tr1 = 1 (%)                     2702 (24.7)   0.359       0.021
  log_dias_treat_imp_sin_na (mean (SD))    4.99 (1.10)   0.074       0.043
  edad_al_ing_1 (mean (SD))               32.82 (8.70)  <0.001       0.542
  ano_nac_corr (mean (SD))              1979.83 (8.99)  <0.001       0.446
  susinidumrec_coc = 1 (%)                  519 ( 4.7)   0.693       0.010
  susinidumrec_pbc = 1 (%)                  872 ( 8.0)  <0.001       0.088
  susinidumrec_mar = 1 (%)                 3865 (35.3)  <0.001       0.341
  susinidumrec_otr = 1 (%)                  238 ( 2.2)   0.569       0.014
  psycom_dum_study = 1 (%)                 2233 (20.4)   0.002       0.071
  psycom_dum_with = 1 (%)                  4850 (44.4)   0.008       0.060
  freq_cons_dum_5day = 1 (%)               5229 (47.8)  <0.001       0.107
  cond_oc_dum_2inact = 1 (%)               1942 (17.8)   0.033       0.048
  cond_oc_dum_3unemp = 1 (%)               4507 (41.2)  <0.001       0.191
  susprindumrec_coc = 1 (%)                2078 (19.0)  <0.001       0.187
  susprindumrec_pbc = 1 (%)                5996 (54.8)  <0.001       0.346
  susprindumrec_mar = 1 (%)                 653 ( 6.0)  <0.001       0.162
  susprindumrec_otr = 1 (%)                 145 ( 1.3)   0.090       0.039
  tipo_de_plan_2_mod (%)                                <0.001       0.298
     basic ambulatory                      3320 (30.4)                    
     GP intensive ambulatory               4212 (38.5)                    
     GP residential                        1906 (17.4)                    
     WO intensive ambulatory                607 ( 5.6)                    
     WO residential                         889 ( 8.1)                    
Code
smd_bin(
  prop.table(table(data_mine_miss_restr_proc2_baseline$susprindum_mar[data_mine_miss_restr_proc2_baseline$policonsumo==0]))[[2]],
  prop.table(table(data_mine_miss_restr_proc2_baseline$susprindum_mar[data_mine_miss_restr_proc2_baseline$policonsumo==1]))[[2]]
)
[1] -0.1630994
Code
smd_bin(
    prop.table(table(data_mine_miss_restr_proc2_baseline$comp_bpsc_y3_severe [data_mine_miss_restr_proc2_baseline$policonsumo==0]))[[2]],
    prop.table(table(data_mine_miss_restr_proc2_baseline$comp_bpsc_y3_severe [data_mine_miss_restr_proc2_baseline$policonsumo==1]))[[2]]
)
[1] -0.3140315
Code
folder_path <- ifelse(dir.exists("E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/"),
                      "E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/",
                      "C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/")

as.data.frame.TableOne(tableOne, smd=T, nonnormal= T)%>% 
    dplyr::mutate(char2=characteristic) %>% 
    tidyr::fill(char2) %>% 
    dplyr::select(char2,everything()) %>% 
    dplyr::mutate(level=ifelse(is.na(level),"[Missing]",level)) %>% 
    dplyr::mutate(char2=dplyr::case_when(characteristic=="NA"~NA_character_,TRUE~as.character(characteristic))) %>% 
    format_cells(1, 1:length(names(.)), "bold") %>%
    dplyr::select(-1) %>% 
    dplyr::mutate_all(~ str_replace_all(., pattern = "\\( ", replacement = "\\(")) %>% 
    dplyr::mutate_all(~ str_trim(.)) %>% 
    dplyr::select(characteristic, level, `0`, `1`, Overall, SMD)%>% 
  {
    knitr::kable(.,size=10, format="html",caption= "Summary descriptives, Polysubstance(1) and no Polysubstance use (0)", escape=T)%>% 
     kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>% 
     kableExtra::scroll_box(width = "100%", height = "375px")  
    write.table(.,paste0(folder_path,"baseline_psu_desc_subsamp_smd.csv"), dec=",", sep="\t")
  }

Time for this code chunk to run: 0 minutes

Code
  rbind.data.frame(
    cbind.data.frame(
      type= "PSU at admission and at least one dropout from the first admission",
    biostat3::survRate(Surv(cens_time, sum_tr_outcomes>0) ~ policonsumo2, 
                       data= data_mine_miss_restr_proc2 %>% 
                        dplyr::group_by(hash_key) %>% 
                        dplyr::mutate(sum_tr_outcomes=sum(tr_outcome==1,na.rm=T)) %>% 
                        dplyr::filter(is_first_occurrence==1) %>% 
                        dplyr::ungroup())  %>% 
      dplyr::mutate(across(c("rate", "lower", "upper"),~sprintf("%1.2f",.*1000)))
    ),
    cbind.data.frame(
      type= "PSU at admission and first dropout",
    biostat3::survRate(Surv(cens_time, tr_outcome==1) ~ policonsumo2, 
                       data= data_mine_miss_restr_proc2 %>% 
                        dplyr::group_by(hash_key) %>% 
                        dplyr::mutate(sum_tr_outcomes=sum(tr_outcome==1,na.rm=T)) %>% 
                        dplyr::filter(is_first_occurrence==1) %>% 
                        dplyr::ungroup())  %>% 
      dplyr::mutate(across(c("rate", "lower", "upper"),~sprintf("%1.2f",.*1000))) 
    ),
    cbind.data.frame(
      type= "At least one treatment reporting PSU and at least one dropout from the first admission",
    biostat3::survRate(Surv(cens_time, sum_tr_outcomes>0) ~ policonsumo2, 
                       data= data_mine_miss_restr_proc2 %>% 
                        dplyr::group_by(hash_key) %>% 
                        dplyr::mutate(sum_tr_outcomes=sum(tr_outcome==1,na.rm=T),
                                      sum_policonsumo2=sum(policonsumo2==1,na.rm=T)) %>% 
                        dplyr::filter(is_first_occurrence==1) %>% 
                        dplyr::ungroup() %>% dplyr::mutate(policonsumo2=ifelse(sum_policonsumo2>0,1,0)))  %>% 
      dplyr::mutate(across(c("rate", "lower", "upper"),~sprintf("%1.2f",.*1000))) 
    ),
    cbind.data.frame(
      type= "At least one treatment reporting PSU and first dropout",
    biostat3::survRate(Surv(cens_time, tr_outcome==1) ~ policonsumo2, 
                       data= data_mine_miss_restr_proc2 %>% 
                        dplyr::group_by(hash_key) %>% 
                        dplyr::mutate(sum_tr_outcomes=sum(tr_outcome==1,na.rm=T),
                                      sum_policonsumo2=sum(policonsumo2==1,na.rm=T)) %>% 
                        dplyr::filter(is_first_occurrence==1) %>% 
                        dplyr::ungroup() %>% dplyr::mutate(policonsumo2=ifelse(sum_policonsumo2>0,1,0)))  %>% 
      dplyr::mutate(across(c("rate", "lower", "upper"),~sprintf("%1.2f",.*1000))) 
    )    
  ) %>%  
    dplyr::mutate(`IR (95% CI)`= paste0(rate, " (", lower,", ", upper,")")) %>% 
    dplyr::select(-any_of(c("rate", "lower", "upper"))) %>% 
    dplyr::mutate(across(c("tstop", "event"),~ format(as.numeric(.), big.mark=","))) %>% 
    
  {
    copiar_nombres2(.)
    write.table(.,paste0(folder_path,"person_time_irrs.csv"), dec=",", sep="\t")
    knitr::kable(.,"html", caption= "Incidence rates (per 1.000 person-months)")%>% 
     kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>% 
     kableExtra::scroll_box(width = "100%", height = "375px")  
  }
Incidence rates (per 1.000 person-months)
type policonsumo2 tstop event IR (95% CI)
policonsumo2=0 PSU at admission and at least one dropout from the first admission 0 161,852.30 2,135 13.19 (12.64, 13.76)
policonsumo2=1 PSU at admission and at least one dropout from the first admission 1 872,863.25 10,085 11.55 (11.33, 11.78)
policonsumo2=01 PSU at admission and first dropout 0 161,852.30 1,833 11.33 (10.81, 11.86)
policonsumo2=11 PSU at admission and first dropout 1 872,863.25 8,615 9.87 (9.66, 10.08)
policonsumo2=02 At least one treatment reporting PSU and at least one dropout from the first admission 0 78,926.85 1,099 13.92 (13.11, 14.77)
policonsumo2=12 At least one treatment reporting PSU and at least one dropout from the first admission 1 955,788.69 11,121 11.64 (11.42, 11.85)
policonsumo2=03 At least one treatment reporting PSU and first dropout 0 78,926.85 936 11.86 (11.11, 12.64)
policonsumo2=13 At least one treatment reporting PSU and first dropout 1 955,788.69 9,512 9.95 (9.75, 10.15)

Time for this code chunk to run: 0 minutes

0.b. IIWs

Code
rbind.data.frame(
  cbind.data.frame(model= "Primary (no corrections, lag=0)", t(matrix(summary(data_mine_miss_restr_proc2$iiw_nocorr_st)))),
  cbind.data.frame(model= "Alternative (no corrections, lag=1)", t(matrix(summary(data_mine_miss_restr_proc2$iiw_nocorr_alt_st)))),
  cbind.data.frame(model= "Primary (after PH, lag=0)", t(matrix(summary(data_mine_miss_restr_proc2_iiw_after_ph$iiw_after_ph_st)))),
  cbind.data.frame(model= "Alternative (after PH, lag=1)", t(matrix(summary(data_mine_miss_restr_proc2_iiw_after_ph_alt$iiw_after_ph_st)))),
  cbind.data.frame(model= "Primary (strata, lag=0)", t(matrix(summary(data_mine_miss_proc2_iiw_strata_alt$iiw_strata_st)))),
  cbind.data.frame(model= "Alternative (strata, lag=1)", t(matrix(summary(data_mine_miss_proc2_iiw_strata_alt_alt$iiw_strata_st))))
) %>%  #Min.    First quartile  Median  Mean    Third quartile  Max.
  dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.)) %>% 
  dplyr::mutate(model= dplyr::case_when(model== "Primary (no corrections, lag=0)"~ "No time-varying coefficients, lagged covariates fixed in 0", model== "Alternative (no corrections, lag=1)"~ "No time-varying coefficients, lagged covariates fixed in 1", model== "Primary (after PH, lag=0)"~ "With time-varying coefficients, lagged covariates fixed in 0", model== "Alternative (after PH, lag=1)"~ "With time-varying coefficients, lagged covariates fixed in 1", model== "Primary (strata, lag=0)"~ "Stratified by follow-up intervals, lagged covariates fixed in 0", model== "Alternative (strata, lag=1)"~ "Stratified by follow-up intervals, lagged covariates fixed in 1")) %>% 
  {
    copiar_nombres2(.)
    write.table(.,paste0(folder_path,"iiws_20240516.csv"), dec=",", sep="\t")
    knitr::kable(.,size=10, format="html",caption= "Descriptive characterization of inverse intensity weights", escape=T, col.names= c("Visit intensity model", "Min.", "First quartile",   "Median",   "Mean", "Third quartile",   "Max."))%>% 
     kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>% 
     kableExtra::scroll_box(width = "100%", height = "375px")  
  }
Descriptive characterization of inverse intensity weights
Visit intensity model Min. First quartile Median Mean Third quartile Max.
No time-varying coefficients, lagged covariates fixed in 0 0.21 1.00 1.62 1.83 2.56 3.79
No time-varying coefficients, lagged covariates fixed in 1 0.13 0.42 0.68 0.70 1.00 1.00
With time-varying coefficients, lagged covariates fixed in 0 0.47 0.67 1.08 1.33 1.90 3.03
With time-varying coefficients, lagged covariates fixed in 1 0.10 0.17 0.55 0.88 1.36 3.56
Stratified by follow-up intervals, lagged covariates fixed in 0 0.29 0.60 0.86 0.99 1.18 3.02
Stratified by follow-up intervals, lagged covariates fixed in 1 0.12 0.12 0.44 0.60 0.76 3.07

Time for this code chunk to run: 0 minutes

0.c. Intensity model

Code
#Cox Counting process


#https://rdrr.io/github/gforge/Greg/src/R/tidy.rms.R
source(paste0(folder_path,"tidy_rms.R"))

term_mapping <- data.frame(
  term = c("lag_tr_outcome", "lag_comp_bpsc_y3_severe", "lag_less_90d_tr1", "log_lag_dias_treat_imp_sin_na", "lag_policonsumo2", "edad_al_ing_1", "ano_nac_corr", "susinidumrec_coc", "susinidumrec_pbc", "susinidumrec_mar", "susinidumrec_otr", "psycom_dum_study", "psycom_dum_with", "freq_cons_dum_5day", "cond_oc_dum_2inact", "cond_oc_dum_3unemp", "susprindumrec_coc", "susprindumrec_pbc", "susprindumrec_mar", "susprindumrec_otr",  "lag_tr_outcome_rec", "lag_less_90d_tr1_rec", "lag_comp_bpsc_y3_severe_rec", "susinidum_coc_rec2", "psycom_dum_with_rec2"),
  term_label = c("Treatment outcome of the previous treatment", "Previous biopsychosocial compromise (severe)", 
                 "Previous treatment duration (<90 days)", "Previous treatment duration (in logarithmic scaled days)", 
                 "Polysubstance use status of the previous treatment", "Age at admission to treatment", "Birth year", 
                  "Primary substance (initial diagnosis), cocaine", "Primary substance (initial diagnosis), cocaine base paste", "Primary substance (initial diagnosis), marijuana", "Primary substance (initial diagnosis), other", "Psychiatric comorbidity (diagnosis unknown or under study)", "Psychiatric comorbidity (confirmed comorbidity)", "Daily frequence of primary substance use at admission", "Occupational status (inactive)", "Occupational status (unemployed)", "Primary substance at admission to treatment (cocaine hydrochloride)", "Primary substance at admission to treatment (cocaine base paste)", "Primary substance at admission to treatment (marijuana)", "Primary substance at admission to treatment (other)", "Treatment outcome of the previous treatment (recoded for interaction with time)", "Previous biopsychosocial compromise (severe, recoded for interaction with time)", "Previous treatment duration (<90 days, recoded for interaction with time)", "Primary substance (initial diagnosis), cocaine (recoded for interaction with time)", "Psychiatric comorbidity (confirmed comorbidity, recoded for interaction with time)")
)

cph_nocorr<-
  cph(Surv(lag_time,time,event)~
        cluster(id)+ 
        lag_tr_outcome+
        lag_comp_bpsc_y3_severe+
        lag_less_90d_tr1+
        log_lag_dias_treat_imp_sin_na +
        lag_policonsumo2 + 
        edad_al_ing_1 + 
        ano_nac_corr + 
        susinidumrec_otr +
        susinidumrec_coc +
        susinidumrec_pbc +
        susinidumrec_mar +
        psycom_dum_study +
        psycom_dum_with +
        freq_cons_dum_5day +
        cond_oc_dum_2inact +
        cond_oc_dum_3unemp +
        susprindumrec_coc +
        susprindumrec_pbc +
        susprindumrec_mar +
        susprindumrec_otr +    
        strat(tipo_de_plan_2_mod),
      data= data_mine_miss_restr_proc2exp %>% data.table::as.data.table() %>% data.frame(), 
      x=TRUE, y=TRUE, surv=TRUE, iter.max = 250*4, tol = 1e-6)

cph_after_ph<-
cph(Surv(lag_time,time,event==1)~ 
                  cluster(id)+ 
                  lag_tr_outcome_rec +
                  log_lag_dias_treat_imp_sin_na +
                  lag_less_90d_tr1_rec+
                  lag_comp_bpsc_y3_severe_rec + 
                  lag_policonsumo2 + 
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidum_coc_rec2 +
                  susinidumrec_otr +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_with_rec2 +
                  psycom_dum_study + 
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                  susprindumrec_coc +
                  susprindumrec_pbc +
                  susprindumrec_mar +
                  susprindumrec_otr +    
                  strat(tipo_de_plan_2_mod), 
                data=data_mine_miss_restr_proc2exp, 
                x=TRUE, y=TRUE, surv=TRUE, iter.max = 250*4, tol = 1e-6)

model_after_time_strat<-
cph(Surv(lag_time,time,event)~
          cluster(id)+ 
          lag_tr_outcome+
          lag_comp_bpsc_y3_severe+
          lag_less_90d_tr1+
          log_lag_dias_treat_imp_sin_na +
          lag_policonsumo2 + 
          edad_al_ing_1 + 
          ano_nac_corr + 
          susinidumrec_otr +
          susinidumrec_coc +
          susinidumrec_pbc +
          susinidumrec_mar +
          psycom_dum_study +
          psycom_dum_with +
          freq_cons_dum_5day +
          cond_oc_dum_2inact +
          cond_oc_dum_3unemp +
          susprindumrec_coc +
          susprindumrec_pbc +
          susprindumrec_mar +
          susprindumrec_otr +
          strat(tipo_de_plan_2_mod)+
          strat(time_interval3),
    data= data_mine_miss_restr_proc2exp %>% data.table::as.data.table() %>% data.frame(), 
    x=TRUE, y=TRUE, surv=TRUE, iter.max = 250*4, tol = 1e-6)

model_after_time_strat_alt<-
cph(Surv(lag_time,time,event)~
          cluster(id)+ 
          lag_tr_outcome+
          lag_comp_bpsc_y3_severe+
          lag_less_90d_tr1+
          log_lag_dias_treat_imp_sin_na +
          lag_policonsumo2 + 
          edad_al_ing_1 + 
          ano_nac_corr + 
          susinidumrec_otr +
          susinidumrec_coc +
          susinidumrec_pbc +
          susinidumrec_mar +
          psycom_dum_study +
          psycom_dum_with +
          freq_cons_dum_5day +
          cond_oc_dum_2inact +
          cond_oc_dum_3unemp +
          susprindumrec_otr +
          susprindumrec_coc +
          susprindumrec_pbc +
          susprindumrec_mar +
          strat(tipo_de_plan_2_mod)+
          strat(time_interval3_alt),
    data= data_mine_miss_restr_proc2exp %>% data.table::as.data.table() %>% data.frame(), 
    x=TRUE, y=TRUE, surv=TRUE, iter.max = 250*4, tol = 1e-6)

invisible("Problemático")
# broom::tidy(coxph(Surv(lag_time,time,event==1)~ 
#         cluster(id)+ 
#         lag_tr_outcome_rec +
#         log_lag_dias_treat_imp_sin_na +
#         lag_less_90d_tr1_rec+
#         lag_comp_bpsc_y3_severe_rec + 
#         lag_policonsumo2 + 
#         edad_al_ing_1 + 
#         ano_nac_corr + 
#         susinidum_coc_rec2 +
#         susinidumrec_otr +
#         susinidum_pbc +
#         susinidum_mar +
#         psycom_dum_with_rec2 +
#         psycom_dum_study + 
#         freq_cons_dum_5day +
#         cond_oc_dum_2inact +
#         cond_oc_dum_3unemp +
#         susprindum_oh +
#         susprindum_coc +
#         susprindum_pbc +
#         susprindum_mar+
#         strata(tipo_de_plan_2_mod), 
#     data=data_mine_miss_restr_proc2), exponentiate=T, conf.int=T)

models_hr_table2<-
rbind.data.frame(
cbind.data.frame(model="No correction", tidy.rms(cph_nocorr, conf.int=T, exponentiate=T)[, c("term","estimate", "conf.low", "conf.high")]%>% 
  dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.))),
cbind.data.frame(model="After PH", tidy.rms(cph_after_ph, conf.int=T, exponentiate=T)[, c("term","estimate", "conf.low", "conf.high")]%>% 
  dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.))),
cbind.data.frame(model="Stratifying", tidy.rms(model_after_time_strat, conf.int=T, exponentiate=T)[, c("term","estimate", "conf.low", "conf.high")]%>% 
  dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.))),
cbind.data.frame(model="Stratifying, alt strata", tidy.rms(model_after_time_strat_alt, conf.int=T, exponentiate=T)[, c("term","estimate", "conf.low", "conf.high")] %>% 
  dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.)))
)

models_hr_table2$term_label <- term_mapping$term_label[match(models_hr_table2$term, term_mapping$term)]

models_hr_table2%>%
  dplyr::select(term_label, everything()) %>% 
  {
    #copiar_nombres2(.)
    write.table(.,paste0(folder_path,"coxph_intensity_model_prev.csv"), dec=",", sep="\t")
    knitr::kable(.,size=10, format="html",caption= "Intensity model", escape=T)%>% 
     kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>% 
     kableExtra::scroll_box(width = "100%", height = "375px")  
  }
Intensity model
term_label model term estimate conf.low conf.high
Treatment outcome of the previous treatment No correction lag_tr_outcome 1.17 1.13 1.21
Previous biopsychosocial compromise (severe) No correction lag_comp_bpsc_y3_severe 1.06 1.03 1.10
Previous treatment duration (<90 days) No correction lag_less_90d_tr1 1.11 1.06 1.16
Previous treatment duration (in logarithmic scaled days) No correction log_lag_dias_treat_imp_sin_na 0.98 0.96 1.00
Polysubstance use status of the previous treatment No correction lag_policonsumo2 0.98 0.95 1.02
Age at admission to treatment No correction edad_al_ing_1 1.26 1.25 1.27
Birth year No correction ano_nac_corr 1.27 1.26 1.28
Primary substance (initial diagnosis), other No correction susinidumrec_otr 0.96 0.87 1.05
Primary substance (initial diagnosis), cocaine No correction susinidumrec_coc 1.07 0.99 1.15
Primary substance (initial diagnosis), cocaine base paste No correction susinidumrec_pbc 1.12 1.06 1.19
Primary substance (initial diagnosis), marijuana No correction susinidumrec_mar 1.12 1.08 1.16
Psychiatric comorbidity (diagnosis unknown or under study) No correction psycom_dum_study 1.03 0.99 1.08
Psychiatric comorbidity (confirmed comorbidity) No correction psycom_dum_with 1.02 0.99 1.05
Daily frequence of primary substance use at admission No correction freq_cons_dum_5day 1.01 0.98 1.05
Occupational status (inactive) No correction cond_oc_dum_2inact 1.06 1.02 1.11
Occupational status (unemployed) No correction cond_oc_dum_3unemp 1.06 1.03 1.10
Primary substance at admission to treatment (cocaine hydrochloride) No correction susprindumrec_coc 0.99 0.95 1.04
Primary substance at admission to treatment (cocaine base paste) No correction susprindumrec_pbc 1.01 0.97 1.05
Primary substance at admission to treatment (marijuana) No correction susprindumrec_mar 0.99 0.92 1.06
Primary substance at admission to treatment (other) No correction susprindumrec_otr 1.11 0.98 1.25
Treatment outcome of the previous treatment (recoded for interaction with time) After PH lag_tr_outcome_rec 0.09 0.08 0.10
Previous treatment duration (in logarithmic scaled days) After PH log_lag_dias_treat_imp_sin_na 0.90 0.87 0.92
Previous biopsychosocial compromise (severe, recoded for interaction with time) After PH lag_less_90d_tr1_rec 0.62 0.49 0.77
Previous treatment duration (<90 days, recoded for interaction with time) After PH lag_comp_bpsc_y3_severe_rec 0.19 0.16 0.22
Polysubstance use status of the previous treatment After PH lag_policonsumo2 0.99 0.95 1.03
Age at admission to treatment After PH edad_al_ing_1 1.20 1.19 1.21
Birth year After PH ano_nac_corr 1.21 1.20 1.22
Primary substance (initial diagnosis), cocaine (recoded for interaction with time) After PH susinidum_coc_rec2 1.04 0.94 1.15
Primary substance (initial diagnosis), other After PH susinidumrec_otr 1.00 0.89 1.12
Primary substance (initial diagnosis), cocaine base paste After PH susinidumrec_pbc 1.14 1.07 1.21
Primary substance (initial diagnosis), marijuana After PH susinidumrec_mar 1.13 1.09 1.17
Psychiatric comorbidity (confirmed comorbidity, recoded for interaction with time) After PH psycom_dum_with_rec2 1.05 1.01 1.09
Psychiatric comorbidity (diagnosis unknown or under study) After PH psycom_dum_study 1.17 1.12 1.23
Daily frequence of primary substance use at admission After PH freq_cons_dum_5day 1.05 1.02 1.09
Occupational status (inactive) After PH cond_oc_dum_2inact 1.08 1.03 1.13
Occupational status (unemployed) After PH cond_oc_dum_3unemp 1.10 1.06 1.14
Primary substance at admission to treatment (cocaine hydrochloride) After PH susprindumrec_coc 0.98 0.93 1.04
Primary substance at admission to treatment (cocaine base paste) After PH susprindumrec_pbc 1.05 1.00 1.10
Primary substance at admission to treatment (marijuana) After PH susprindumrec_mar 0.95 0.88 1.03
Primary substance at admission to treatment (other) After PH susprindumrec_otr 1.00 0.87 1.14
Treatment outcome of the previous treatment Stratifying lag_tr_outcome 1.15 1.11 1.19
Previous biopsychosocial compromise (severe) Stratifying lag_comp_bpsc_y3_severe 1.05 1.02 1.09
Previous treatment duration (<90 days) Stratifying lag_less_90d_tr1 1.11 1.06 1.17
Previous treatment duration (in logarithmic scaled days) Stratifying log_lag_dias_treat_imp_sin_na 0.98 0.96 1.00
Polysubstance use status of the previous treatment Stratifying lag_policonsumo2 0.98 0.94 1.02
Age at admission to treatment Stratifying edad_al_ing_1 1.22 1.21 1.23
Birth year Stratifying ano_nac_corr 1.22 1.21 1.23
Primary substance (initial diagnosis), other Stratifying susinidumrec_otr 0.96 0.87 1.05
Primary substance (initial diagnosis), cocaine Stratifying susinidumrec_coc 1.07 0.99 1.15
Primary substance (initial diagnosis), cocaine base paste Stratifying susinidumrec_pbc 1.10 1.04 1.17
Primary substance (initial diagnosis), marijuana Stratifying susinidumrec_mar 1.12 1.09 1.16
Psychiatric comorbidity (diagnosis unknown or under study) Stratifying psycom_dum_study 1.01 0.97 1.06
Psychiatric comorbidity (confirmed comorbidity) Stratifying psycom_dum_with 1.01 0.98 1.05
Daily frequence of primary substance use at admission Stratifying freq_cons_dum_5day 1.01 0.98 1.04
Occupational status (inactive) Stratifying cond_oc_dum_2inact 1.08 1.03 1.13
Occupational status (unemployed) Stratifying cond_oc_dum_3unemp 1.09 1.05 1.13
Primary substance at admission to treatment (cocaine hydrochloride) Stratifying susprindumrec_coc 1.00 0.96 1.05
Primary substance at admission to treatment (cocaine base paste) Stratifying susprindumrec_pbc 1.02 0.98 1.06
Primary substance at admission to treatment (marijuana) Stratifying susprindumrec_mar 1.00 0.93 1.08
Primary substance at admission to treatment (other) Stratifying susprindumrec_otr 1.08 0.95 1.24
Treatment outcome of the previous treatment Stratifying, alt strata lag_tr_outcome 1.05 1.01 1.09
Previous biopsychosocial compromise (severe) Stratifying, alt strata lag_comp_bpsc_y3_severe 1.02 0.99 1.06
Previous treatment duration (<90 days) Stratifying, alt strata lag_less_90d_tr1 1.04 0.99 1.09
Previous treatment duration (in logarithmic scaled days) Stratifying, alt strata log_lag_dias_treat_imp_sin_na 0.97 0.96 0.99
Polysubstance use status of the previous treatment Stratifying, alt strata lag_policonsumo2 0.99 0.95 1.03
Age at admission to treatment Stratifying, alt strata edad_al_ing_1 1.05 1.05 1.06
Birth year Stratifying, alt strata ano_nac_corr 1.05 1.05 1.06
Primary substance (initial diagnosis), other Stratifying, alt strata susinidumrec_otr 1.03 0.93 1.14
Primary substance (initial diagnosis), cocaine Stratifying, alt strata susinidumrec_coc 1.10 1.02 1.18
Primary substance (initial diagnosis), cocaine base paste Stratifying, alt strata susinidumrec_pbc 1.02 0.97 1.08
Primary substance (initial diagnosis), marijuana Stratifying, alt strata susinidumrec_mar 1.02 0.99 1.06
Psychiatric comorbidity (diagnosis unknown or under study) Stratifying, alt strata psycom_dum_study 1.00 0.96 1.05
Psychiatric comorbidity (confirmed comorbidity) Stratifying, alt strata psycom_dum_with 0.98 0.95 1.01
Daily frequence of primary substance use at admission Stratifying, alt strata freq_cons_dum_5day 1.02 0.99 1.05
Occupational status (inactive) Stratifying, alt strata cond_oc_dum_2inact 0.99 0.94 1.03
Occupational status (unemployed) Stratifying, alt strata cond_oc_dum_3unemp 1.02 0.99 1.06
Primary substance at admission to treatment (other) Stratifying, alt strata susprindumrec_otr 0.92 0.80 1.06
Primary substance at admission to treatment (cocaine hydrochloride) Stratifying, alt strata susprindumrec_coc 0.98 0.93 1.03
Primary substance at admission to treatment (cocaine base paste) Stratifying, alt strata susprindumrec_pbc 0.98 0.94 1.03
Primary substance at admission to treatment (marijuana) Stratifying, alt strata susprindumrec_mar 0.99 0.92 1.06

Time for this code chunk to run: 0.1 minutes

0.c.Intensity model (with model that originated IWWs)

Taking other substances as dummy just as in irrelong-21-iiw-weights-nocorr (step 3).

Code
term_mapping2 <- data.frame(
  term = c("lag_tr_outcome", "lag_comp_bpsc_y3_severe", "lag_less_90d_tr1", "log_lag_dias_treat_imp_sin_na", "lag_policonsumo2", "edad_al_ing_1", "ano_nac_corr", "susinidum_coc", "susinidum_pbc", "susinidum_mar", "susinidum_oh", "psycom_dum_study", "psycom_dum_with", "freq_cons_dum_5day", "cond_oc_dum_2inact", "cond_oc_dum_3unemp", "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh",  "lag_tr_outcome_rec", "lag_less_90d_tr1_rec", "lag_comp_bpsc_y3_severe_rec", "susinidum_coc_rec2", "psycom_dum_with_rec2"),
  term_label = c("Treatment outcome of the previous treatment", "Previous biopsychosocial compromise (severe)", 
                 "Previous treatment duration (<90 days)", "Previous treatment duration (in logarithmic scaled days)", 
                 "Polysubstance use status of the previous treatment", "Age at admission to treatment", "Birth year", 
                  "Primary substance (initial diagnosis), cocaine", "Primary substance (initial diagnosis), cocaine base paste", "Primary substance (initial diagnosis), marijuana", "Primary substance (initial diagnosis), alcohol", "Psychiatric comorbidity (diagnosis unknown or under study)", "Psychiatric comorbidity (confirmed comorbidity)", "Daily frequence of primary substance use at admission", "Occupational status (inactive)", "Occupational status (unemployed)", "Primary substance at admission to treatment (cocaine hydrochloride)", "Primary substance at admission to treatment (cocaine base paste)", "Primary substance at admission to treatment (marijuana)", "Primary substance at admission to treatment (alcohol)", "Treatment outcome of the previous treatment (recoded for interaction with time)", "Previous biopsychosocial compromise (severe, recoded for interaction with time)", "Previous treatment duration (<90 days, recoded for interaction with time)", "Primary substance (initial diagnosis), cocaine (recoded for interaction with time)", "Psychiatric comorbidity (confirmed comorbidity, recoded for interaction with time)")
)

#susinidumrec_oh susprindumrec_oh

cph_nocorr_oh<-
  cph(Surv(lag_time,time,event)~
        cluster(id)+ 
        lag_tr_outcome+ #tr_outcome.lag
        lag_comp_bpsc_y3_severe+ #comp_bpsc_y3_severe.lag
        lag_less_90d_tr1+ #less_90d_tr1.lag
        log_lag_dias_treat_imp_sin_na + #log_dias_treat_imp_sin_na.lag
        lag_policonsumo2 +  #policonsumo2.lag 
        edad_al_ing_1 + 
        ano_nac_corr + 
        susinidum_oh+
        #susinidumrec_otr +
        susinidum_coc +
        susinidum_pbc +
        susinidum_mar +
        psycom_dum_study +
        psycom_dum_with +
        freq_cons_dum_5day +
        cond_oc_dum_2inact +
        cond_oc_dum_3unemp +
        susprindum_oh +
        susprindum_coc +
        susprindum_pbc +
        susprindum_mar +
        strat(tipo_de_plan_2_mod),
      data= data_mine_miss_restr_proc2 %>% data.table::as.data.table() %>% data.frame(), 
      x=TRUE, y=TRUE, surv=TRUE, iter.max = 250*4, tol = 1e-6)

models_hr_table2_oh<-
rbind.data.frame(
cbind.data.frame(model="No correction", tidy.rms(cph_nocorr_oh, conf.int=T, exponentiate=T)[, c("term","estimate", "conf.low", "conf.high")]%>% dplyr::mutate_if(is.numeric, ~sprintf("%1.2f",.)))
)

models_hr_table2_oh$term_label <- term_mapping2$term_label[match(models_hr_table2_oh$term, term_mapping2$term)]

models_hr_table2_oh%>%
  dplyr::select(term_label, everything()) %>% 
  {
    #copiar_nombres2(.)
    write.table(.,paste0(folder_path,"coxph_intensity_model_prev_250216.csv"), dec=",", sep="\t")
    knitr::kable(.,size=10, format="html",caption= "Intensity model (no correction, original)", escape=T)%>% 
     kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>% 
     kableExtra::scroll_box(width = "100%", height = "375px")  
  }
Intensity model (no correction, original)
term_label model term estimate conf.low conf.high
Treatment outcome of the previous treatment No correction lag_tr_outcome 1.17 1.13 1.21
Previous biopsychosocial compromise (severe) No correction lag_comp_bpsc_y3_severe 1.06 1.03 1.10
Previous treatment duration (<90 days) No correction lag_less_90d_tr1 1.11 1.06 1.16
Previous treatment duration (in logarithmic scaled days) No correction log_lag_dias_treat_imp_sin_na 0.98 0.96 1.00
Polysubstance use status of the previous treatment No correction lag_policonsumo2 0.98 0.95 1.02
Age at admission to treatment No correction edad_al_ing_1 1.26 1.25 1.27
Birth year No correction ano_nac_corr 1.27 1.26 1.28
Primary substance (initial diagnosis), alcohol No correction susinidum_oh 1.05 0.95 1.14
Primary substance (initial diagnosis), cocaine No correction susinidum_coc 1.12 1.00 1.25
Primary substance (initial diagnosis), cocaine base paste No correction susinidum_pbc 1.17 1.06 1.30
Primary substance (initial diagnosis), marijuana No correction susinidum_mar 1.17 1.07 1.29
Psychiatric comorbidity (diagnosis unknown or under study) No correction psycom_dum_study 1.03 0.99 1.08
Psychiatric comorbidity (confirmed comorbidity) No correction psycom_dum_with 1.02 0.99 1.05
Daily frequence of primary substance use at admission No correction freq_cons_dum_5day 1.01 0.98 1.05
Occupational status (inactive) No correction cond_oc_dum_2inact 1.06 1.02 1.11
Occupational status (unemployed) No correction cond_oc_dum_3unemp 1.06 1.03 1.10
Primary substance at admission to treatment (alcohol) No correction susprindum_oh 0.90 0.80 1.02
Primary substance at admission to treatment (cocaine hydrochloride) No correction susprindum_coc 0.90 0.79 1.01
Primary substance at admission to treatment (cocaine base paste) No correction susprindum_pbc 0.91 0.81 1.03
Primary substance at admission to treatment (marijuana) No correction susprindum_mar 0.89 0.78 1.02

Time for this code chunk to run: 0 minutes

0.d. Cox intensity model behind weights, under two scenarios (lag1 and lag0)

Code
update_terms <- function(tidy_table, term_mapping2) {
  
  # Function to convert tidy_table term to mapping key.
  convert_term <- function(term) {
    # Check if the term ends with ".lag"
    if(grepl("\\.lag$", term)) {
      base <- sub("\\.lag$", "", term)
      # Special case: if the base is "log_dias_treat_imp_sin_na", insert "lag" after "log_"
      if(base == "log_dias_treat_imp_sin_na") {
        return("log_lag_dias_treat_imp_sin_na")
      } else {
        return(paste0("lag_", base))
      }
    } else {
      return(term)
    }
  }
  
  # Create a new column with the mapping key.
  tidy_table$term_key <- sapply(tidy_table$term, convert_term)
  
  # Merge the tidy_table with the term_mapping on the mapping key.
  # term_mapping has columns "term" and "term_label".
  tidy_table <- merge(tidy_table, term_mapping2, by.x = "term_key", by.y = "term", all.x = TRUE)
  
  # Replace the original term with the label if available.
  tidy_table$term <- ifelse(!is.na(tidy_table$term_label), tidy_table$term_label, tidy_table$term)
  
  # Remove temporary columns used for the merge.
  tidy_table$term_key <- NULL
  tidy_table$term_label <- NULL
  
  return(tidy_table)
}

tidy_table <- tidy_coxph(iiw_nocorr$m)

tidy_table_alt <- tidy_coxph(iiw_nocorr_alt$m)

term_vector <- c(
  "Treatment outcome of the previous treatment",
  "Previous treatment duration (in logarithmic scaled days)",
  "Previous treatment duration (<90 days)",
  "Previous biopsychosocial compromise (severe)",
  "Polysubstance use status of the previous treatment",
  "Age at admission to treatment",
  "Birth year",
  "Primary substance (initial diagnosis), cocaine",
  "Primary substance (initial diagnosis), alcohol",
  "Primary substance (initial diagnosis), cocaine base paste",
  "Primary substance (initial diagnosis), marijuana",
  "Psychiatric comorbidity (confirmed comorbidity)",
  "Psychiatric comorbidity (diagnosis unknown or under study)",
  "Daily frequence of primary substance use at admission",
  "Occupational status (inactive)",
  "Occupational status (unemployed)",
  "Primary substance at admission to treatment (alcohol)",
  "Primary substance at admission to treatment (cocaine hydrochloride)",
  "Primary substance at admission to treatment (cocaine base paste)",
  "Primary substance at admission to treatment (marijuana)"
)

rbind.data.frame(cbind.data.frame(lag="Lag 0", update_terms(format_hr_ci(tidy_table, digits = 2), term_mapping2)),cbind.data.frame(lag="Lag 1", update_terms(format_hr_ci(tidy_table_alt, digits = 2), term_mapping2)))%>% 
  dplyr::mutate(term= factor(term, levels=term_vector))%>% 
  dplyr::arrange(lag, term)%>% #rio::export("clipboard")
  {
    #copiar_nombres2(.)
    write.table(., paste0(folder_path,"coxph_intensity_model_250216.csv"), dec=",", sep="\t")
    knitr::kable(tidyr::pivot_wider(dplyr::select(.,"term","lag", "hr_ci"), names_from="lag", values_from="hr_ci"), size=10, format="html",caption= "Specifications of the treatment (visit) process, in Hazard Ratios (HR) and under different lagged variables scenarios (no correction, original)", escape=T)%>% 
     kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>% 
     kableExtra::scroll_box(width = "100%", height = "375px")  
  }
Specifications of the treatment (visit) process, in Hazard Ratios (HR) and under different lagged variables scenarios (no correction, original)
term Lag 0 Lag 1
Treatment outcome of the previous treatment 0.65 (0.63, 0.67) 1.30 (1.26, 1.34)
Previous treatment duration (in logarithmic scaled days) 0.85 (0.84, 0.86) 1.48 (1.44, 1.53)
Previous treatment duration (<90 days) 0.69 (0.67, 0.72) 2.63 (2.52, 2.74)
Previous biopsychosocial compromise (severe) 0.88 (0.85, 0.90) 1.54 (1.50, 1.57)
Polysubstance use status of the previous treatment 0.62 (0.61, 0.64) 1.26 (1.22, 1.30)
Age at admission to treatment 1.05 (1.05, 1.06) 1.06 (1.06, 1.06)
Birth year 1.06 (1.06, 1.06) 1.06 (1.06, 1.06)
Primary substance (initial diagnosis), cocaine 1.01 (0.95, 1.08) 1.05 (0.98, 1.12)
Primary substance (initial diagnosis), alcohol 1.01 (0.96, 1.07) 1.04 (0.99, 1.10)
Primary substance (initial diagnosis), cocaine base paste 0.98 (0.92, 1.04) 1.05 (0.99, 1.11)
Primary substance (initial diagnosis), marijuana 1.04 (0.99, 1.10) 1.03 (0.98, 1.09)
Psychiatric comorbidity (confirmed comorbidity) 1.00 (0.99, 1.02) 0.95 (0.94, 0.97)
Psychiatric comorbidity (diagnosis unknown or under study) 1.00 (0.98, 1.03) 0.74 (0.72, 0.76)
Daily frequence of primary substance use at admission 1.02 (1.01, 1.04) 0.95 (0.94, 0.97)
Occupational status (inactive) 0.99 (0.96, 1.01) 0.96 (0.94, 0.99)
Occupational status (unemployed) 1.03 (1.01, 1.05) 0.98 (0.96, 1.00)
Primary substance at admission to treatment (alcohol) 0.99 (0.93, 1.06) 1.03 (0.96, 1.09)
Primary substance at admission to treatment (cocaine hydrochloride) 1.11 (1.04, 1.19) 1.03 (0.97, 1.10)
Primary substance at admission to treatment (cocaine base paste) 1.14 (1.07, 1.22) 1.01 (0.95, 1.08)
Primary substance at admission to treatment (marijuana) 1.04 (0.96, 1.12) 1.05 (0.98, 1.12)

Time for this code chunk to run: 0 minutes

GEE

1. No correction

Code
plan_names <- attr(table(data_mine_miss_restr_proc2$tipo_de_plan_2_mod),"dimnames")[[1]]

Time for this code chunk to run: 0 minutes

1.1. No correciton, no weight

Initialize a list to store models

Code
nocorr_nowgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geese(tr_outcome ~ policonsumo2 +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, #origen_ingreso_mod fis_comorbidity_icd_10 dg_trs_cons_sus_or
                 id = id, 
                 data = current_data,
                 family = poisson(), 
                 corstr = "independence", 
                 jack = T)
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  nocorr_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}

Time for this code chunk to run: 0 minutes

1.1.1 No correction, no weight, GEEM Poisson

Initialize a list to store models

Code
geem_nocorr_nowgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset using geem
  model <- geem(tr_outcome ~ policonsumo2 +
                  comp_bpsc_y3_severe+
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidumrec_otr +
                  susinidumrec_coc +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_study +
                  psycom_dum_with +
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                  susprindumrec_coc +
                  susprindumrec_pbc +
                  susprindumrec_mar+
                  susprindumrec_otr, # Include relevant predictors
                id = id, 
                data = current_data,
                family = poisson("log"), 
                corstr = "independence")
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  if(model$converged==FALSE){stop("model did not coverge")}
  cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
  geem_nocorr_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}
QIC for plan basic ambulatory = 19453.6
QIC for plan GP intensive ambulatory = 22261.9
QIC for plan GP residential = 9818.1
QIC for plan WO intensive ambulatory = 3424.2
QIC for plan WO residential = 4826.1

Time for this code chunk to run: 0 minutes

1.1.2 No correction, no weight, GEEM Negative binomial

Define a sequence of x values

Code
x_values <- seq(0.2, 16, by = 0.1)

Time for this code chunk to run: 0 minutes

Initialize a list to store models and their QICs

Code
geem2_nocorr_nowgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit models

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
  
  # Initialize variable to store the best QIC and associated model
  best_qic <- Inf
  best_model <- NULL
  best_x <- NA
  
  # Loop through each x value to find the one with the lowest QIC
  for (x in x_values) {
    # Fit the GEE model for the current subset using geem with negative binomial
    model <- geem(tr_outcome ~ policonsumo2 +
                    comp_bpsc_y3_severe+
                    edad_al_ing_1 + 
                    ano_nac_corr + 
                    susinidumrec_otr +
                    susinidumrec_coc +
                    susinidumrec_pbc +
                    susinidumrec_mar +
                    psycom_dum_study +
                    psycom_dum_with +
                    freq_cons_dum_5day +
                    cond_oc_dum_2inact +
                    cond_oc_dum_3unemp +
                    susprindumrec_coc +
                    susprindumrec_pbc +
                    susprindumrec_mar+
                    susprindumrec_otr, # Include relevant predictors
                  id = id, 
                  data = current_data,
                  family = MASS::negative.binomial(x),
                  corstr = "independence")
    
    if (model$converged) {
      # Calculate QIC for the model
      model_qic <- MuMIn::QIC(model)
      
      # Check if this model has a better QIC
      if (model_qic < best_qic) {
        best_qic <- model_qic
        best_model <- model
        best_x <- x
      }
    }
  }
  
  # Store the best model and its details in the list
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  geem2_nocorr_nowgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
  
  # Optional: Print out some information on progress
  cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}
Best QIC for plan basic ambulatory with x = 16.0: 19838.564776
Best QIC for plan GP intensive ambulatory with x = 16.0: 22686.829620
Best QIC for plan GP residential with x = 16.0: 9979.192804
Best QIC for plan WO intensive ambulatory with x = 16.0: 3487.845045
Best QIC for plan WO residential with x = 16.0: 4909.067252

Time for this code chunk to run: 0 minutes

1.2. No correction, weight

Initialize a list to store models

Code
nocorr_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geese(tr_outcome ~ policonsumo2 +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                 id = id, 
                 data = current_data,
                 family = poisson(), 
                 weight = current_data$iiw_nocorr_st,
                 corstr = "independence", 
                 jack = T)
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  nocorr_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}

Time for this code chunk to run: 0 minutes

1.2.1 No correction, no weight, GEEM Poisson

Initialize a list to store models

Code
geem_nocorr_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset using geem
  model <- geem(tr_outcome ~ policonsumo2 +
                  comp_bpsc_y3_severe+
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidumrec_otr +
                  susinidumrec_coc +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_study +
                  psycom_dum_with +
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                  susprindumrec_coc +
                  susprindumrec_pbc +
                  susprindumrec_mar +
                  susprindumrec_otr, # Include relevant predictors
                id = id, 
                data = current_data,
                weight = current_data$iiw_nocorr_st,
                family = poisson("log"), 
                corstr = "independence")
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  if(model$converged==FALSE){stop("model did not coverge")}
  cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
  geem_nocorr_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}
QIC for plan basic ambulatory = 19465.1
QIC for plan GP intensive ambulatory = 22275.9
QIC for plan GP residential = 9833.7
QIC for plan WO intensive ambulatory = 3434.7
QIC for plan WO residential = 4839.1

Time for this code chunk to run: 0 minutes

1.2.2 No correction, no weight, GEEM Negative binomial

Define a sequence of x values

Code
x_values <- seq(158e4, 159e4, by = 100)

Time for this code chunk to run: 0 minutes

Initialize a list to store models and their QICs

Code
geem2_nocorr_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit models

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
  
  # Initialize variable to store the best QIC and associated model
  best_qic <- Inf
  best_model <- NULL
  best_x <- NA
  
  # Loop through each x value to find the one with the lowest QIC
  for (x in x_values) {
    # Fit the GEE model for the current subset using geem with negative binomial
    model <- geem(tr_outcome ~ policonsumo2 +
                    comp_bpsc_y3_severe+
                    edad_al_ing_1 + 
                    ano_nac_corr + 
                    susinidumrec_otr +
                    susinidumrec_coc +
                    susinidumrec_pbc +
                    susinidumrec_mar +
                    psycom_dum_study +
                    psycom_dum_with +
                    freq_cons_dum_5day +
                    cond_oc_dum_2inact +
                    cond_oc_dum_3unemp +
                    susprindumrec_coc +
                    susprindumrec_pbc +
                    susprindumrec_mar +
                    susprindumrec_otr, # Include relevant predictors
                  id = id, 
                  data = current_data,
                  weight = current_data$iiw_nocorr_st,
                  family = MASS::negative.binomial(x),
                  corstr = "independence")
    
    if (model$converged) {
      # Calculate QIC for the model
      model_qic <- MuMIn::QIC(model)
      
      # Check if this model has a better QIC
      if (model_qic < best_qic) {
        best_qic <- model_qic
        best_model <- model
        best_x <- x
      }
    }
  }
  
  # Store the best model and its details in the list
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  geem2_nocorr_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
  
  # Optional: Print out some information on progress
  cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}
Best QIC for plan basic ambulatory with x = 1588300.0: 19465.056444
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22275.954182
Best QIC for plan GP residential with x = 1588300.0: 9833.688227
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3434.680800
Best QIC for plan WO residential with x = 1588300.0: 4839.094573

Time for this code chunk to run: 0 minutes

Check the result for a particular plan

Code
print(geem2_nocorr_wgt_list[[1]])
$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe + 
    edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc + 
    susinidumrec_pbc + susinidumrec_mar + psycom_dum_study + 
    psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact + 
    cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc + 
    susprindumrec_mar + susprindumrec_otr, id = id, data = current_data, 
    family = MASS::negative.binomial(x), corstr = "independence", 
    weights = current_data$iiw_nocorr_st)

 Coefficients: 
 (Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
       -9.22      0.02195            -0.03892     -0.002341     0.004535
 susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
          0.03033        -0.004384          0.09993          0.02421
 psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
           0.1066         0.03428             0.0419           -0.05004
 cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
           0.005096           0.01384           0.04043          0.002462
 susprindumrec_otr
           -0.1842

 Scale Parameter:  0.4019 

 Correlation Model:  independence
 Estimated Correlation Parameter:  0 

 Number of clusters:  4360   Maximum cluster size:  10 
 Number of observations with nonzero weight:  9993 

$best_x
[1] 1588300

$best_qic
     QIC 
19465.06 

Time for this code chunk to run: 0 minutes

1.3. No correction, alt. weight

Initialize a list to store models

Code
nocorr_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geese(tr_outcome ~ policonsumo2 +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                 id = id, 
                 data = current_data,
                 family = poisson(), 
                 weight = current_data$iiw_nocorr_alt_st,
                 corstr = "independence", 
                 jack = T)
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  nocorr_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}

Time for this code chunk to run: 0 minutes

1.3.1 No correction, alt weight, GEEM Poisson

Initialize a list to store models

Code
geem_nocorr_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset using geem
  model <- geem(tr_outcome ~ policonsumo2 +
                  comp_bpsc_y3_severe+
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidumrec_otr +
                  susinidumrec_coc +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_study +
                  psycom_dum_with +
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                  susprindumrec_coc +
                  susprindumrec_pbc +
                  susprindumrec_mar +
                  susprindumrec_otr, # Include relevant predictors
                id = id, 
                data = current_data,
                weight = current_data$iiw_nocorr_alt_st,
                family = poisson("log"), 
                corstr = "independence")
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  if(model$converged==FALSE){stop("model did not coverge")}
  cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
  geem_nocorr_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}
QIC for plan basic ambulatory = 19458.8
QIC for plan GP intensive ambulatory = 22271.9
QIC for plan GP residential = 9827.8
QIC for plan WO intensive ambulatory = 3431.3
QIC for plan WO residential = 4835.6

Time for this code chunk to run: 0 minutes

1.3.2 No correction, alt weight, GEEM Negative binomial

Code
x_values <- seq(158e4, 159e4, by = 100)

Time for this code chunk to run: 0 minutes

Initialize a list to store models and their QICs

Code
geem2_nocorr_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit models

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp, tipo_de_plan_2_mod == plan_names[i])
  
  # Initialize variable to store the best QIC and associated model
  best_qic <- Inf
  best_model <- NULL
  best_x <- NA
  
  # Loop through each x value to find the one with the lowest QIC
  for (x in x_values) {
    # Fit the GEE model for the current subset using geem with negative binomial
    model <- geem(tr_outcome ~ policonsumo2 +
                    comp_bpsc_y3_severe+
                    edad_al_ing_1 + 
                    ano_nac_corr + 
                    susinidumrec_otr +
                    susinidumrec_coc +
                    susinidumrec_pbc +
                    susinidumrec_mar +
                    psycom_dum_study +
                    psycom_dum_with +
                    freq_cons_dum_5day +
                    cond_oc_dum_2inact +
                    cond_oc_dum_3unemp +
                    susprindumrec_coc +
                    susprindumrec_pbc +
                    susprindumrec_mar +
                    susprindumrec_otr, # Include relevant predictors
                  id = id, 
                  data = current_data,
                  weight = current_data$iiw_nocorr_alt_st,
                  family = MASS::negative.binomial(x),
                  corstr = "independence")
    
    if (model$converged) {
      # Calculate QIC for the model
      model_qic <- MuMIn::QIC(model)
      
      # Check if this model has a better QIC
      if (model_qic < best_qic) {
        best_qic <- model_qic
        best_model <- model
        best_x <- x
      }
    }
  }
  
  # Store the best model and its details in the list
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  geem2_nocorr_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
  
  # Optional: Print out some information on progress
  cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}
Best QIC for plan basic ambulatory with x = 1588300.0: 19458.771969
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22271.862242
Best QIC for plan GP residential with x = 1588300.0: 9827.753103
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3431.276563
Best QIC for plan WO residential with x = 1588300.0: 4835.575200

Time for this code chunk to run: 0 minutes

Check the result for a particular plan

Code
print(geem2_nocorr_alt_wgt_list[[1]])
$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe + 
    edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc + 
    susinidumrec_pbc + susinidumrec_mar + psycom_dum_study + 
    psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact + 
    cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc + 
    susprindumrec_mar + susprindumrec_otr, id = id, data = current_data, 
    family = MASS::negative.binomial(x), corstr = "independence", 
    weights = current_data$iiw_nocorr_alt_st)

 Coefficients: 
 (Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
       -19.1      0.02437            -0.03728      0.003479     0.009432
 susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
          0.07601        -0.005822           0.0957          0.02467
 psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
           0.1475         0.04159            0.03779           -0.04085
 cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
            0.01445           0.01352           0.04788          0.007057
 susprindumrec_otr
           -0.1853

 Scale Parameter:  0.1465 

 Correlation Model:  independence
 Estimated Correlation Parameter:  0 

 Number of clusters:  4360   Maximum cluster size:  10 
 Number of observations with nonzero weight:  9993 

$best_x
[1] 1588300

$best_qic
     QIC 
19458.77 

Time for this code chunk to run: 0 minutes

2. After PH

2.1. After PH, no weight

Initialize a list to store models

Code
afterph_nowgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
data_mine_miss_restr_proc2exp_iiw_after_ph <-
  data_mine_miss_restr_proc2_iiw_after_ph %>% 
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>% 
  {if (nrow(.) > nrow(data_mine_miss_restr_proc2_iiw_after_ph)) {message("left join added rows"); .} else .} %>% 
  dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>% 
  dplyr::mutate(
    susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
    susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
    susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
    susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
    susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
    susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
    susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
    susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
  )

for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geese(tr_outcome ~ policonsumo2 +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                 id = id, 
                 data = current_data,
                 family = poisson(), 
                 #weight= current_data$iiw_nocorr_alt_st,
                 corstr = "independence", 
                 jack = T)
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  afterph_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}

Time for this code chunk to run: 0 minutes

2.1.1 After PH, no weight, GEEM Poisson

Initialize a list to store models

Code
geem_afterph_nowgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset using geem
  model <- geem(tr_outcome ~ policonsumo2 +
                  comp_bpsc_y3_severe+
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidumrec_otr +
                  susinidumrec_coc +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_study +
                  psycom_dum_with +
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                  susprindumrec_coc +
                  susprindumrec_pbc +
                  susprindumrec_mar +
                  susprindumrec_otr, # Include relevant predictors
                id = id, 
                data = current_data,
                family = poisson("log"), 
                corstr = "independence")
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  if(model$converged==FALSE){stop("model did not coverge")}
  cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
  geem_afterph_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}
QIC for plan basic ambulatory = 19453.6
QIC for plan GP intensive ambulatory = 22261.9
QIC for plan GP residential = 9818.1
QIC for plan WO intensive ambulatory = 3424.2
QIC for plan WO residential = 4826.1

Time for this code chunk to run: 0 minutes

2.1.2 After PH, no weight, GEEM Negative binomial

Define a sequence of x values

Code
x_values <- seq(158e4, 159e4, by = 100)

Time for this code chunk to run: 0 minutes

Initialize a list to store models and their QICs

Code
geem2_afterph_nowgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit models

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
  
  # Initialize variable to store the best QIC and associated model
  best_qic <- Inf
  best_model <- NULL
  best_x <- NA
  
  # Loop through each x value to find the one with the lowest QIC
  for (x in x_values) {
    # Fit the GEE model for the current subset using geem with negative binomial
    model <- geem(tr_outcome ~ policonsumo2 +
                    comp_bpsc_y3_severe+
                    edad_al_ing_1 + 
                    ano_nac_corr + 
                    susinidumrec_otr +
                    susinidumrec_coc +
                    susinidumrec_pbc +
                    susinidumrec_mar +
                    psycom_dum_study +
                    psycom_dum_with +
                    freq_cons_dum_5day +
                    cond_oc_dum_2inact +
                    cond_oc_dum_3unemp +
                    susprindumrec_coc +
                    susprindumrec_pbc +
                    susprindumrec_mar +
                    susprindumrec_otr, # Include relevant predictors
                  id = id, 
                  data = current_data,
                  family = MASS::negative.binomial(x),
                  corstr = "independence")
    
    if (model$converged) {
      # Calculate QIC for the model
      model_qic <- MuMIn::QIC(model)
      
      # Check if this model has a better QIC
      if (model_qic < best_qic) {
        best_qic <- model_qic
        best_model <- model
        best_x <- x
      }
    }
  }
  
  # Store the best model and its details in the list
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  geem2_afterph_nowgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
  
  # Optional: Print out some information on progress
  cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}
Best QIC for plan basic ambulatory with x = 1588300.0: 19453.628013
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22261.947187
Best QIC for plan GP residential with x = 1588300.0: 9818.145865
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3424.239894
Best QIC for plan WO residential with x = 1588300.0: 4826.145839

Time for this code chunk to run: 0 minutes

Check the result for a particular plan

Code
print(geem2_nocorr_nowgt_list[[1]])
$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe + 
    edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc + 
    susinidumrec_pbc + susinidumrec_mar + psycom_dum_study + 
    psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact + 
    cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc + 
    susprindumrec_mar + susprindumrec_otr, id = id, data = current_data, 
    family = MASS::negative.binomial(x), corstr = "independence")

 Coefficients: 
 (Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
      -12.32      0.02507            -0.03872     -0.000283     0.006061
 susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
          0.05421        -0.003447          0.09881          0.02315
 psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
           0.1383         0.04429            0.04057           -0.04427
 cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
             0.0111           0.01587           0.04783          0.007784
 susprindumrec_otr
           -0.1904

 Scale Parameter:  0.2033 

 Correlation Model:  independence
 Estimated Correlation Parameter:  0 

 Number of clusters:  4360   Maximum cluster size:  10 
 Number of observations with nonzero weight:  9993 

$best_x
[1] 16

$best_qic
     QIC 
19838.56 

Time for this code chunk to run: 0 minutes

As expected, QIC is equal, because the covariates and observations were the same

2.2. After PH, weight

Initialize a list to store models

Code
afterph_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geese(tr_outcome ~ policonsumo2 +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                 id = id, 
                 data = current_data,
                 family = poisson(), 
                 weight = current_data$iiw_after_ph_st,
                 corstr = "independence", 
                 jack = T)
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  afterph_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}

Time for this code chunk to run: 0 minutes

2.2.1 After PH, weight, GEEM Poisson

Initialize a list to store models

Code
geem_afterph_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset using geem
  model <- geem(tr_outcome ~ policonsumo2 +
                  comp_bpsc_y3_severe+
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidumrec_otr +
                  susinidumrec_coc +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_study +
                  psycom_dum_with +
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                  susprindumrec_coc +
                  susprindumrec_pbc +
                  susprindum_mar +
                  susprindumrec_otr, # Include relevant predictors
                id = id, 
                data = current_data,
                weight = current_data$iiw_after_ph_st,
                family = poisson("log"), 
                corstr = "independence")
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  if(model$converged==FALSE){stop("model did not coverge")}
  cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
  geem_afterph_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}
QIC for plan basic ambulatory = 19479.4
QIC for plan GP intensive ambulatory = 22284.4
QIC for plan GP residential = 9830.3
QIC for plan WO intensive ambulatory = 3445.9
QIC for plan WO residential = 4839.7

Time for this code chunk to run: 0 minutes

2.2.2 After PH, weight, GEEM Negative binomial

Define a sequence of x values

Code
x_values <- seq(158e4, 159e4, by = 100)

Time for this code chunk to run: 0 minutes

Initialize a list to store models and their QICs

Code
geem2_afterph_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit models

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph, tipo_de_plan_2_mod == plan_names[i])
  
  # Initialize variable to store the best QIC and associated model
  best_qic <- Inf
  best_model <- NULL
  best_x <- NA
  
  # Loop through each x value to find the one with the lowest QIC
  for (x in x_values) {
    # Fit the GEE model for the current subset using geem with negative binomial
    model <- geem(tr_outcome ~ policonsumo2 +
                    comp_bpsc_y3_severe+
                    edad_al_ing_1 + 
                    ano_nac_corr + 
                    susinidumrec_otr +
                    susinidumrec_coc +
                    susinidumrec_pbc +
                    susinidumrec_mar +
                    psycom_dum_study +
                    psycom_dum_with +
                    freq_cons_dum_5day +
                    cond_oc_dum_2inact +
                    cond_oc_dum_3unemp +
                    susprindumrec_coc +
                    susprindumrec_pbc +
                    susprindumrec_mar +
                    susprindumrec_otr, # Include relevant predictors
                  id = id, 
                  data = current_data,
                  weight = current_data$iiw_after_ph_st,
                  family = MASS::negative.binomial(x),
                  corstr = "independence")
    
    if (model$converged) {
      # Calculate QIC for the model
      model_qic <- MuMIn::QIC(model)
      
      # Check if this model has a better QIC
      if (model_qic < best_qic) {
        best_qic <- model_qic
        best_model <- model
        best_x <- x
      }
    }
  }
  
  # Store the best model and its details in the list
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  geem2_afterph_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
  
  # Optional: Print out some information on progress
  cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}
Best QIC for plan basic ambulatory with x = 1588300.0: 19479.434224
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22284.436679
Best QIC for plan GP residential with x = 1588300.0: 9830.254882
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3445.890741
Best QIC for plan WO residential with x = 1588300.0: 4839.713265

Time for this code chunk to run: 0 minutes

Check the result for a particular plan

Code
print(geem2_afterph_wgt_list[[1]])
$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe + 
    edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc + 
    susinidumrec_pbc + susinidumrec_mar + psycom_dum_study + 
    psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact + 
    cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc + 
    susprindumrec_mar + susprindumrec_otr, id = id, data = current_data, 
    family = MASS::negative.binomial(x), corstr = "independence", 
    weights = current_data$iiw_after_ph_st)

 Coefficients: 
 (Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
      -5.745      0.02667            -0.04339     -0.004904      0.00281
 susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
           0.0327         -0.01607           0.1156          0.02389
 psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
           0.1134          0.0387            0.04731           -0.05169
 cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
           0.002699           0.01144           0.03865          0.003677
 susprindumrec_otr
           -0.1873

 Scale Parameter:  0.2933 

 Correlation Model:  independence
 Estimated Correlation Parameter:  0 

 Number of clusters:  4360   Maximum cluster size:  10 
 Number of observations with nonzero weight:  9993 

$best_x
[1] 1588300

$best_qic
     QIC 
19479.43 

Time for this code chunk to run: 0 minutes

2.3. After PH, alt. weight

Initialize a list to store models

Code
afterph_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
data_mine_miss_restr_proc2exp_iiw_after_ph_alt <-
  data_mine_miss_restr_proc2_iiw_after_ph_alt %>% 
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>% 
  {if (nrow(.) > nrow(data_mine_miss_restr_proc2_iiw_after_ph_alt)) {message("left join added rows"); .} else .} %>% 
  dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>% 
  dplyr::mutate(
    susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
    susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
    susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
    susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
    susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
    susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
    susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
    susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
  )


for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph_alt, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geese(tr_outcome ~ policonsumo2 +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                 id = id, 
                 data = current_data,
                 family = poisson(), 
                 weight = current_data$iiw_after_ph_st,
                 corstr = "independence", 
                 jack = T)
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  afterph_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}

Time for this code chunk to run: 0 minutes

2.3.1 After PH, alt weight, GEEM Poisson

Initialize a list to store models

Code
geem_afterph_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph_alt, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset using geem
  model <- geem(tr_outcome ~ policonsumo2 +
                  comp_bpsc_y3_severe+
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidumrec_otr +
                  susinidumrec_coc +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_study +
                  psycom_dum_with +
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                  susprindumrec_coc +
                  susprindumrec_pbc +
                  susprindumrec_mar +
                  susprindumrec_otr, # Include relevant predictors
                id = id, 
                data = current_data,
                weight = current_data$iiw_after_ph_st,
                family = poisson("log"), 
                corstr = "independence")
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  if(model$converged==FALSE){stop("model did not coverge")}
  cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
  geem_afterph_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}
QIC for plan basic ambulatory = 19519.4
QIC for plan GP intensive ambulatory = 22324.6
QIC for plan GP residential = 9868.4
QIC for plan WO intensive ambulatory = 3468.6
QIC for plan WO residential = 4861.0

Time for this code chunk to run: 0 minutes

2.3.2 After PH, alt weight, GEEM Negative binomial

Define a sequence of x values

Code
x_values <- seq(158e4, 159e4, by = 100)

Time for this code chunk to run: 0 minutes

Initialize a list to store models and their QICs

Code
geem2_afterph_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit models

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp_iiw_after_ph_alt, tipo_de_plan_2_mod == plan_names[i])
  
  # Initialize variable to store the best QIC and associated model
  best_qic <- Inf
  best_model <- NULL
  best_x <- NA
  
  # Loop through each x value to find the one with the lowest QIC
  for (x in x_values) {
    # Fit the GEE model for the current subset using geem with negative binomial
    model <- geem(tr_outcome ~ policonsumo2 +
                    comp_bpsc_y3_severe+
                    edad_al_ing_1 + 
                    ano_nac_corr + 
                    susinidumrec_otr +
                    susinidumrec_coc +
                    susinidumrec_pbc +
                    susinidumrec_mar +
                    psycom_dum_study +
                    psycom_dum_with +
                    freq_cons_dum_5day +
                    cond_oc_dum_2inact +
                    cond_oc_dum_3unemp +
                    susprindumrec_coc +
                    susprindumrec_pbc +
                    susprindumrec_mar +
                    susprindumrec_otr, # Include relevant predictors
                  id = id, 
                  data = current_data,
                  weight = current_data$iiw_after_ph_st,
                  family = MASS::negative.binomial(x),
                  corstr = "independence")
    
    if (model$converged) {
      # Calculate QIC for the model
      model_qic <- MuMIn::QIC(model)
      
      # Check if this model has a better QIC
      if (model_qic < best_qic) {
        best_qic <- model_qic
        best_model <- model
        best_x <- x
      }
    }
  }
  
  # Store the best model and its details in the list
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  geem2_afterph_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
  
  # Optional: Print out some information on progress
  cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}
Best QIC for plan basic ambulatory with x = 1588300.0: 19519.425246
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22324.624851
Best QIC for plan GP residential with x = 1588300.0: 9868.351832
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3468.579298
Best QIC for plan WO residential with x = 1588300.0: 4861.024120

Time for this code chunk to run: 0 minutes

Check the result for a particular plan

Code
print(geem2_afterph_alt_wgt_list[[1]])
$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe + 
    edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc + 
    susinidumrec_pbc + susinidumrec_mar + psycom_dum_study + 
    psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact + 
    cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc + 
    susprindumrec_mar + susprindumrec_otr, id = id, data = current_data, 
    family = MASS::negative.binomial(x), corstr = "independence", 
    weights = current_data$iiw_after_ph_st)

 Coefficients: 
 (Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
      0.7589      0.02652            -0.03843     -0.007929   -0.0004416
 susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
          0.03431         -0.04666           0.1411          0.03227
 psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
          0.09837         0.04042             0.0413           -0.06383
 cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
           0.004229           0.03314            0.0559           0.03698
 susprindumrec_otr
           -0.1908

 Scale Parameter:  0.2291 

 Correlation Model:  independence
 Estimated Correlation Parameter:  0 

 Number of clusters:  4360   Maximum cluster size:  10 
 Number of observations with nonzero weight:  9993 

$best_x
[1] 1588300

$best_qic
     QIC 
19519.43 

Time for this code chunk to run: 0 minutes

3. Stratifying

3.1. Stratifying, no weight

Initialize a list to store models

Code
strata_nowgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
data_mine_miss_proc2exp_iiw_strata <-
  data_mine_miss_proc2_iiw_strata %>% 
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>% 
  {if (nrow(.) > nrow(data_mine_miss_restr_proc2_iiw_after_ph_alt)) {message("left join added rows"); .} else .} %>% 
  dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>% 
  dplyr::mutate(
    susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
    susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
    susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
    susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
    susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
    susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
    susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
    susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
  )

for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geese(tr_outcome ~ policonsumo2 +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                 id = id, 
                 data = current_data,
                 family = poisson(), 
                 #weight= current_data$iiw_nocorr_alt_st,
                 corstr = "independence", 
                 jack = T)
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  strata_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}

Time for this code chunk to run: 0 minutes

3.1.1 Stratifying, no weight, GEEM Poisson

Initialize a list to store models

Code
geem_strata_nowgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset using geem
  model <- geem(tr_outcome ~ policonsumo2 +
                  comp_bpsc_y3_severe+
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidumrec_otr +
                  susinidumrec_coc +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_study +
                  psycom_dum_with +
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                  susprindumrec_coc +
                  susprindumrec_pbc +
                  susprindumrec_mar +
                  susprindumrec_otr, # Include relevant predictors
                id = id, 
                data = current_data,
                family = poisson("log"), 
                corstr = "independence")
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  if(model$converged==FALSE){stop("model did not coverge")}
  cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
  geem_strata_nowgt_list[[paste("model", model_name, sep = "_")]] <- model
}
QIC for plan basic ambulatory = 19453.6
QIC for plan GP intensive ambulatory = 22261.9
QIC for plan GP residential = 9818.1
QIC for plan WO intensive ambulatory = 3424.2
QIC for plan WO residential = 4826.1

Time for this code chunk to run: 0 minutes

3.1.2 Stratifying, no weight, GEEM Negative binomial

Define a sequence of x values

Code
x_values <- seq(158e4, 159e4, by = 100)

Time for this code chunk to run: 0 minutes

Initialize a list to store models and their QICs

Code
geem2_strata_nowgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit models

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
  
  # Initialize variable to store the best QIC and associated model
  best_qic <- Inf
  best_model <- NULL
  best_x <- NA
  
  # Loop through each x value to find the one with the lowest QIC
  for (x in x_values) {
    # Fit the GEE model for the current subset using geem with negative binomial
    model <- geem(tr_outcome ~ policonsumo2 +
                    comp_bpsc_y3_severe+
                    edad_al_ing_1 + 
                    ano_nac_corr + 
                    susinidumrec_otr +
                    susinidumrec_coc +
                    susinidumrec_pbc +
                    susinidumrec_mar +
                    psycom_dum_study +
                    psycom_dum_with +
                    freq_cons_dum_5day +
                    cond_oc_dum_2inact +
                    cond_oc_dum_3unemp +
                    susprindumrec_coc +
                    susprindumrec_pbc +
                    susprindumrec_mar +
                    susprindumrec_otr, # Include relevant predictors
                  id = id, 
                  data = current_data,
                  family = MASS::negative.binomial(x),
                  corstr = "independence")
    
    if (model$converged) {
      # Calculate QIC for the model
      model_qic <- MuMIn::QIC(model)
      
      # Check if this model has a better QIC
      if (model_qic < best_qic) {
        best_qic <- model_qic
        best_model <- model
        best_x <- x
      }
    }
  }
  
  # Store the best model and its details in the list
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  geem2_strata_nowgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
  
  # Optional: Print out some information on progress
  cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}
Best QIC for plan basic ambulatory with x = 1588300.0: 19453.628013
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22261.947187
Best QIC for plan GP residential with x = 1588300.0: 9818.145865
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3424.239894
Best QIC for plan WO residential with x = 1588300.0: 4826.145839

Time for this code chunk to run: 0 minutes

Check the result for a particular plan

Code
print(geem2_nocorr_nowgt_list[[1]])
$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe + 
    edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc + 
    susinidumrec_pbc + susinidumrec_mar + psycom_dum_study + 
    psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact + 
    cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc + 
    susprindumrec_mar + susprindumrec_otr, id = id, data = current_data, 
    family = MASS::negative.binomial(x), corstr = "independence")

 Coefficients: 
 (Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
      -12.32      0.02507            -0.03872     -0.000283     0.006061
 susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
          0.05421        -0.003447          0.09881          0.02315
 psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
           0.1383         0.04429            0.04057           -0.04427
 cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
             0.0111           0.01587           0.04783          0.007784
 susprindumrec_otr
           -0.1904

 Scale Parameter:  0.2033 

 Correlation Model:  independence
 Estimated Correlation Parameter:  0 

 Number of clusters:  4360   Maximum cluster size:  10 
 Number of observations with nonzero weight:  9993 

$best_x
[1] 16

$best_qic
     QIC 
19838.56 

Time for this code chunk to run: 0 minutes

3.2. Stratifying, weight

Initialize a list to store models

Code
strata_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geese(tr_outcome ~ policonsumo2 +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                 id = id, 
                 data = current_data,
                 family = poisson(), 
                 weight = current_data$iiw_strata_st,
                 corstr = "independence", 
                 jack = T)
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  strata_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}

Time for this code chunk to run: 0 minutes

3.2.1 Stratifying, weight, GEEM Poisson

Initialize a list to store models

Code
geem_strata_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset using geem
  model <- geem(tr_outcome ~ policonsumo2 +
                  comp_bpsc_y3_severe+
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidumrec_otr +
                  susinidumrec_coc +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_study +
                  psycom_dum_with +
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                id = id, 
                data = current_data,
                weight = current_data$iiw_strata_st,
                family = poisson("log"), 
                corstr = "independence")
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  if(model$converged==FALSE){stop("model did not coverge")}
  cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
  geem_strata_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}
QIC for plan basic ambulatory = 19476.8
QIC for plan GP intensive ambulatory = 22283.0
QIC for plan GP residential = 9839.4
QIC for plan WO intensive ambulatory = 3453.0
QIC for plan WO residential = 4851.1

Time for this code chunk to run: 0 minutes

3.2.2 Stratifying, weight, GEEM Negative binomial

Define a sequence of x values

Code
x_values <- seq(158e4, 159e4, by = 100)

Time for this code chunk to run: 0 minutes

Initialize a list to store models and their QICs

Code
geem2_strata_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit models

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata, tipo_de_plan_2_mod == plan_names[i])
  
  # Initialize variable to store the best QIC and associated model
  best_qic <- Inf
  best_model <- NULL
  best_x <- NA
  
  # Loop through each x value to find the one with the lowest QIC
  for (x in x_values) {
    # Fit the GEE model for the current subset using geem with negative binomial
    model <- geem(tr_outcome ~ policonsumo2 +
                    comp_bpsc_y3_severe+
                    edad_al_ing_1 + 
                    ano_nac_corr + 
                    susinidumrec_otr +
                    susinidumrec_coc +
                    susinidumrec_pbc +
                    susinidumrec_mar +
                    psycom_dum_study +
                    psycom_dum_with +
                    freq_cons_dum_5day +
                    cond_oc_dum_2inact +
                    cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                  id = id, 
                  data = current_data,
                  weight= current_data$iiw_strata_st,
                  family = MASS::negative.binomial(x),
                  corstr = "independence")
    
    if (model$converged) {
      # Calculate QIC for the model
      model_qic <- MuMIn::QIC(model)
      
      # Check if this model has a better QIC
      if (model_qic < best_qic) {
        best_qic <- model_qic
        best_model <- model
        best_x <- x
      }
    }
  }
  
  # Store the best model and its details in the list
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  geem2_strata_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
  
  # Optional: Print out some information on progress
  cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}
Best QIC for plan basic ambulatory with x = 1588300.0: 19476.842070
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22283.006125
Best QIC for plan GP residential with x = 1588300.0: 9839.424427
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3452.975453
Best QIC for plan WO residential with x = 1588300.0: 4851.130363

Time for this code chunk to run: 0 minutes

Check the result for a particular plan

Code
print(geem2_strata_wgt_list[[1]])
$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe + 
    edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc + 
    susinidumrec_pbc + susinidumrec_mar + psycom_dum_study + 
    psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact + 
    cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc + 
    susprindumrec_mar + susprindumrec_otr, id = id, data = current_data, 
    family = MASS::negative.binomial(x), corstr = "independence", 
    weights = current_data$iiw_strata_st)

 Coefficients: 
 (Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
      -15.18      0.01948            -0.05564     0.0002688       0.0075
 susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
          0.05049         0.005513           0.1008          0.02269
 psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
           0.1429         0.05033            0.05562           -0.05253
 cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
            0.01156         -0.001358           0.04337           0.04282
 susprindumrec_otr
           -0.1754

 Scale Parameter:  0.3845 

 Correlation Model:  independence
 Estimated Correlation Parameter:  0 

 Number of clusters:  4360   Maximum cluster size:  10 
 Number of observations with nonzero weight:  9993 

$best_x
[1] 1588300

$best_qic
     QIC 
19476.84 

Time for this code chunk to run: 0 minutes

3.3. Stratifying, alt. weight

Initialize a list to store models

Code
strata_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
data_mine_miss_proc2exp_iiw_strata_alt <-
  data_mine_miss_proc2_iiw_strata_alt %>% 
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>% 
  {if (nrow(.) > nrow(data_mine_miss_restr_proc2_iiw_after_ph_alt)) {message("left join added rows"); .} else .} %>% 
  dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>% 
  dplyr::mutate(
    susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
    susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
    susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
    susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
    susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
    susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
    susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
    susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
  )


for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geese(tr_outcome ~ policonsumo2 +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                 id = id, 
                 data = current_data,
                 family = poisson(), 
                 weight = current_data$iiw_strata_st,
                 corstr = "independence", 
                 jack = T)
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  strata_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}

Time for this code chunk to run: 0 minutes

3.3.1 Stratifying, alt weight, GEEM Poisson

Initialize a list to store models

Code
geem_strata_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset using geem
  model <- geem(tr_outcome ~ policonsumo2 +
                  comp_bpsc_y3_severe+
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidumrec_otr +
                  susinidumrec_coc +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_study +
                  psycom_dum_with +
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                id = id, 
                data = current_data,
                weight = current_data$iiw_strata_st,
                family = poisson("log"), 
                corstr = "independence")
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  if(model$converged==FALSE){stop("model did not coverge")}
  cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
  geem_strata_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}
QIC for plan basic ambulatory = 19467.2
QIC for plan GP intensive ambulatory = 22269.4
QIC for plan GP residential = 9837.5
QIC for plan WO intensive ambulatory = 3437.3
QIC for plan WO residential = 4841.8
Code
#Alternative, better than main

Time for this code chunk to run: 0 minutes

Code
print(geem_strata_alt_wgt_list[[1]])
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe + 
    edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc + 
    susinidumrec_pbc + susinidumrec_mar + psycom_dum_study + 
    psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact + 
    cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc + 
    susprindumrec_mar + susprindumrec_otr, id = id, data = current_data, 
    family = poisson("log"), corstr = "independence", weights = current_data$iiw_strata_st)

 Coefficients: 
 (Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
       -9.62      0.03745             -0.0414     -0.001121     0.004702
 susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
          0.09628         -0.01782           0.1065          0.02583
 psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
           0.1531         0.05717            0.03588           -0.03803
 cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
            0.01523           0.02253           0.04334           0.01237
 susprindumrec_otr
           -0.2359

 Scale Parameter:  0.1917 

 Correlation Model:  independence
 Estimated Correlation Parameter:  0 

 Number of clusters:  4360   Maximum cluster size:  10 
 Number of observations with nonzero weight:  9993 

Time for this code chunk to run: 0 minutes

3.3.2 No correction, alt weight, GEEM Negative binomial

Define a sequence of x values

Code
x_values <- seq(158e4, 159e4, by = 100)

Time for this code chunk to run: 0 minutes

Initialize a list to store models and their QICs

Code
geem2_strata_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit models

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt, tipo_de_plan_2_mod == plan_names[i])
  
  # Initialize variable to store the best QIC and associated model
  best_qic <- Inf
  best_model <- NULL
  best_x <- NA
  
  # Loop through each x value to find the one with the lowest QIC
  for (x in x_values) {
    # Fit the GEE model for the current subset using geem with negative binomial
    model <- geem(tr_outcome ~ policonsumo2 +
                    comp_bpsc_y3_severe+
                    edad_al_ing_1 + 
                    ano_nac_corr + 
                    susinidumrec_otr +
                    susinidumrec_coc +
                    susinidumrec_pbc +
                    susinidumrec_mar +
                    psycom_dum_study +
                    psycom_dum_with +
                    freq_cons_dum_5day +
                    cond_oc_dum_2inact +
                    cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                  id = id, 
                  data = current_data,
                  weight = current_data$iiw_strata_st,
                  family = MASS::negative.binomial(x),
                  corstr = "independence")
    
    if (model$converged) {
      # Calculate QIC for the model
      model_qic <- MuMIn::QIC(model)
      
      # Check if this model has a better QIC
      if (model_qic < best_qic) {
        best_qic <- model_qic
        best_model <- model
        best_x <- x
      }
    }
  }
  
  # Store the best model and its details in the list
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  geem2_strata_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
  
  # Optional: Print out some information on progress
  cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}
Best QIC for plan basic ambulatory with x = 1588300.0: 19467.162599
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22269.414488
Best QIC for plan GP residential with x = 1588300.0: 9837.476338
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3437.325576
Best QIC for plan WO residential with x = 1588300.0: 4841.777059

Time for this code chunk to run: 0 minutes

Check the result for a particular plan

Code
print(geem2_strata_alt_wgt_list[[1]])
$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe + 
    edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc + 
    susinidumrec_pbc + susinidumrec_mar + psycom_dum_study + 
    psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact + 
    cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc + 
    susprindumrec_mar + susprindumrec_otr, id = id, data = current_data, 
    family = MASS::negative.binomial(x), corstr = "independence", 
    weights = current_data$iiw_strata_st)

 Coefficients: 
 (Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
       -9.62      0.03745             -0.0414     -0.001121     0.004702
 susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
          0.09628         -0.01782           0.1065          0.02583
 psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
           0.1531         0.05717            0.03588           -0.03803
 cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
            0.01523           0.02253           0.04334           0.01237
 susprindumrec_otr
           -0.2359

 Scale Parameter:  0.1917 

 Correlation Model:  independence
 Estimated Correlation Parameter:  0 

 Number of clusters:  4360   Maximum cluster size:  10 
 Number of observations with nonzero weight:  9993 

$best_x
[1] 1588300

$best_qic
     QIC 
19467.16 

Time for this code chunk to run: 0 minutes

3.4. Stratifying, 2nd alt. weight

Initialize a list to store models

Code
strata_alt_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
data_mine_miss_proc2exp_iiw_strata_alt_alt <-
  data_mine_miss_proc2_iiw_strata_alt_alt %>% 
left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>% 
  {if (nrow(.) > nrow(data_mine_miss_restr_proc2_iiw_after_ph_alt)) {message("left join added rows"); .} else .} %>% 
  dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>% 
  dplyr::mutate(
    susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
    susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
    susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
    susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
    susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
    susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
    susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
    susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
  )


for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt_alt, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geese(tr_outcome ~ policonsumo2 +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                 id = id, 
                 data = current_data,
                 family = poisson(), 
                 weight = current_data$iiw_strata_st,
                 corstr = "independence", 
                 jack = T)
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  strata_alt_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}

Time for this code chunk to run: 0 minutes

3.4.1 Stratifying, alt weight, GEEM Poisson

Initialize a list to store models

Code
geem_strata_alt_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit a model

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt_alt, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset using geem
  model <- geem(tr_outcome ~ policonsumo2 +
                  comp_bpsc_y3_severe+
                  edad_al_ing_1 + 
                  ano_nac_corr + 
                  susinidumrec_otr +
                  susinidumrec_coc +
                  susinidumrec_pbc +
                  susinidumrec_mar +
                  psycom_dum_study +
                  psycom_dum_with +
                  freq_cons_dum_5day +
                  cond_oc_dum_2inact +
                  cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                id = id, 
                data = current_data,
                weight = current_data$iiw_strata_st,
                family = poisson("log"), 
                corstr = "independence")
  
  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  if(model$converged==FALSE){stop("model did not coverge")}
  cat(sprintf("QIC for plan %s = %.1f\n", plan_names[i], as.numeric(MuMIn::QIC(model))))
  geem_strata_alt_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- model
}
QIC for plan basic ambulatory = 19505.4
QIC for plan GP intensive ambulatory = 22280.0
QIC for plan GP residential = 9885.3
QIC for plan WO intensive ambulatory = 3466.2
QIC for plan WO residential = 4874.2

Time for this code chunk to run: 0 minutes

Code
invisible("Alternative, worst model of strata")

Time for this code chunk to run: 0 minutes

3.4.2 Stratifying, alt weight, GEEM Negative binomial

Define a sequence of x values

Code
x_values <- seq(158e4, 159e4, by = 100)

Time for this code chunk to run: 0 minutes

Initialize a list to store models and their QICs

Code
geem2_strata_alt_alt_wgt_list <- list()

Time for this code chunk to run: 0 minutes

Loop over each unique plan name to fit models

Code
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_proc2exp_iiw_strata_alt_alt, tipo_de_plan_2_mod == plan_names[i])
  
  # Initialize variable to store the best QIC and associated model
  best_qic <- Inf
  best_model <- NULL
  best_x <- NA
  
  # Loop through each x value to find the one with the lowest QIC
  for (x in x_values) {
    # Fit the GEE model for the current subset using geem with negative binomial
    model <- geem(tr_outcome ~ policonsumo2 +
                    comp_bpsc_y3_severe+
                    edad_al_ing_1 + 
                    ano_nac_corr + 
                    susinidumrec_otr +
                    susinidumrec_coc +
                    susinidumrec_pbc +
                    susinidumrec_mar +
                    psycom_dum_study +
                    psycom_dum_with +
                    freq_cons_dum_5day +
                    cond_oc_dum_2inact +
                    cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, 
                  id = id, 
                  data = current_data,
                  weight = current_data$iiw_strata_st,
                  family = MASS::negative.binomial(x),
                  corstr = "independence")
    
    if (model$converged) {
      # Calculate QIC for the model
      model_qic <- MuMIn::QIC(model)
      
      # Check if this model has a better QIC
      if (model_qic < best_qic) {
        best_qic <- model_qic
        best_model <- model
        best_x <- x
      }
    }
  }
  
  # Store the best model and its details in the list
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  geem2_strata_alt_alt_wgt_list[[paste("model", model_name, sep = "_")]] <- list(best_model = best_model, best_x = best_x, best_qic = best_qic)
  
  # Optional: Print out some information on progress
  cat(sprintf("Best QIC for plan %s with x = %.1f: %f\n", plan_names[i], best_x, best_qic))
}
Best QIC for plan basic ambulatory with x = 1588300.0: 19505.358265
Best QIC for plan GP intensive ambulatory with x = 1588300.0: 22279.986625
Best QIC for plan GP residential with x = 1588300.0: 9885.267664
Best QIC for plan WO intensive ambulatory with x = 1588300.0: 3466.220923
Best QIC for plan WO residential with x = 1588300.0: 4874.236495

Time for this code chunk to run: 0 minutes

Check the result for a particular plan

Code
print(geem2_strata_alt_alt_wgt_list[[1]])
$best_model
geem(formula = tr_outcome ~ policonsumo2 + comp_bpsc_y3_severe + 
    edad_al_ing_1 + ano_nac_corr + susinidumrec_otr + susinidumrec_coc + 
    susinidumrec_pbc + susinidumrec_mar + psycom_dum_study + 
    psycom_dum_with + freq_cons_dum_5day + cond_oc_dum_2inact + 
    cond_oc_dum_3unemp + susprindumrec_coc + susprindumrec_pbc + 
    susprindumrec_mar + susprindumrec_otr, id = id, data = current_data, 
    family = MASS::negative.binomial(x), corstr = "independence", 
    weights = current_data$iiw_strata_st)

 Coefficients: 
 (Intercept) policonsumo2 comp_bpsc_y3_severe edad_al_ing_1 ano_nac_corr
      -12.83      0.02767            -0.05236      0.001664      0.00628
 susinidumrec_otr susinidumrec_coc susinidumrec_pbc susinidumrec_mar
           0.1033          -0.0514           0.1021          0.01249
 psycom_dum_study psycom_dum_with freq_cons_dum_5day cond_oc_dum_2inact
           0.1768         0.06588            0.02143           -0.03881
 cond_oc_dum_3unemp susprindumrec_coc susprindumrec_pbc susprindumrec_mar
            0.02918           0.01468           0.04859          0.005722
 susprindumrec_otr
           -0.2477

 Scale Parameter:  0.1097 

 Correlation Model:  independence
 Estimated Correlation Parameter:  0 

 Number of clusters:  4360   Maximum cluster size:  10 
 Number of observations with nonzero weight:  9993 

$best_x
[1] 1588300

$best_qic
     QIC 
19505.36 

Time for this code chunk to run: 0 minutes

Model output

Code
types_of_model<-
  c("No time-varying coefficients, no weight",
    "No time-varying coefficients, lag=0",
    "No time-varying coefficients, lag=1",
    "With time-varying coefficients, no weight",
    "With time-varying coefficients, lag=0",
    "With time-varying coefficients, lag=1",
    "Stratified by follow-up intervals, no weight",
    "Stratified by follow-up intervals, lag=0",
    "Stratified by follow-up intervals, lag=1",
    "Stratified by follow-up intervals, 2nd, lag=0"
)
types_of_model2<-
  c("No time-varying coefficients, no weight",
    "No time-varying coefficients, no weight, NB",
    "No time-varying coefficients, lag=0",
    "No time-varying coefficients, lag=0, NB",
    "No time-varying coefficients, lag=1",
    "No time-varying coefficients, lag=1, NB",
    "With time-varying coefficients, no weight",
    "With time-varying coefficients, no weight, NB",
    "With time-varying coefficients, lag=0",
    "With time-varying coefficients, lag=0, NB",
    "With time-varying coefficients, lag=1",
    "With time-varying coefficients, lag=1, NB",
    "Stratified by follow-up intervals, no weight",
    "Stratified by follow-up intervals, no weight, NB",
    "Stratified by follow-up intervals, lag=0",
    "Stratified by follow-up intervals, lag=0, NB",
    "Stratified by follow-up intervals, lag=1",
    "Stratified by follow-up intervals, lag=1, NB",
    "Stratified by follow-up intervals, 2nd, lag=1",
    "Stratified by follow-up intervals, 2nd, lag=1, NB"    
  )

plan_names[1]
cbind.data.frame(model= types_of_model,
rbind.data.frame(
dplyr::filter(tidy_geese_model(nocorr_nowgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(nocorr_alt_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_nowgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(afterph_alt_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_nowgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_wgt_list[[1]]), Term=="policonsumo2"),
dplyr::filter(tidy_geese_model(strata_alt_alt_wgt_list[[1]]), Term=="policonsumo2")))%>% 
  {
    copiar_nombres2(.)
    assign(paste0("result_data_geese_", gsub(" ","_",plan_names[1])), ., envir = .GlobalEnv)
    print(.)
  }

cbind.data.frame(model= types_of_model2,
     QIC= c(as.numeric(QIC(geem_nocorr_nowgt_list[[1]])),
            as.numeric(geem2_nocorr_nowgt_list[[1]]$best_qic),
            as.numeric(QIC(geem_nocorr_wgt_list[[1]])),
            as.numeric(geem2_nocorr_wgt_list[[1]]$best_qic),
            as.numeric(QIC(geem_nocorr_alt_wgt_list[[1]])),
            as.numeric(geem2_nocorr_alt_wgt_list[[1]]$best_qic),
            as.numeric(QIC(geem_afterph_nowgt_list[[1]])),
            as.numeric(geem2_afterph_nowgt_list[[1]]$best_qic),
            as.numeric(QIC(geem_afterph_wgt_list[[1]])),
            as.numeric(geem2_afterph_wgt_list[[1]]$best_qic),
            as.numeric(QIC(geem_afterph_alt_wgt_list[[1]])),
            as.numeric(geem2_afterph_alt_wgt_list[[1]]$best_qic),
            
            as.numeric(QIC(geem_strata_nowgt_list[[1]])),
            as.numeric(geem2_strata_nowgt_list[[1]]$best_qic),
            as.numeric(QIC(geem_strata_wgt_list[[1]])),
            as.numeric(geem2_strata_wgt_list[[1]]$best_qic),
            as.numeric(QIC(geem_strata_alt_wgt_list[[1]])),
            as.numeric(geem2_strata_alt_wgt_list[[1]]$best_qic),
            as.numeric(QIC(geem_strata_alt_alt_wgt_list[[1]])),
            as.numeric(geem2_strata_alt_alt_wgt_list[[1]]$best_qic)                        
            ),
     rbind.data.frame(
       dplyr::filter(summary2.geem(geem_nocorr_nowgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem2_nocorr_nowgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem_nocorr_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem2_nocorr_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem_nocorr_alt_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem2_nocorr_alt_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem_afterph_nowgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem2_afterph_nowgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem_afterph_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem2_afterph_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem_afterph_alt_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem2_afterph_alt_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
       
       dplyr::filter(summary2.geem(geem_strata_nowgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem2_strata_nowgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),                   
       dplyr::filter(summary2.geem(geem_strata_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem2_strata_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem_strata_alt_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem2_strata_alt_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem_strata_alt_alt_wgt_list[[1]],exponentiate=T, digits=2),Term=="policonsumo2"),
       dplyr::filter(summary2.geem(geem2_strata_alt_alt_wgt_list[[1]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2")
       ))%>% 
  {
    copiar_nombres2(.)
    assign(paste0("result_data_geem_", gsub(" ","_",plan_names[1])), ., envir = .GlobalEnv)
    print(.)
  }

i=2
plan_names[i]
cbind.data.frame(model= types_of_model,
rbind.data.frame(
  dplyr::filter(tidy_geese_model(nocorr_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(nocorr_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(nocorr_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_alt_alt_wgt_list[[i]]), Term=="policonsumo2")))%>% 
  {
    copiar_nombres2(.)  
    print(.)
  }

#[1] "GP intensive ambulatory". significativo
cbind.data.frame(model= types_of_model2,
         QIC= c(as.numeric(QIC(geem_nocorr_nowgt_list[[i]])),
                as.numeric(geem2_nocorr_nowgt_list[[i]]$best_qic),
                as.numeric(QIC(geem_nocorr_wgt_list[[i]])),
                as.numeric(geem2_nocorr_wgt_list[[i]]$best_qic),
                as.numeric(QIC(geem_nocorr_alt_wgt_list[[i]])),
                as.numeric(geem2_nocorr_alt_wgt_list[[i]]$best_qic),
                as.numeric(QIC(geem_afterph_nowgt_list[[i]])),
                as.numeric(geem2_afterph_nowgt_list[[i]]$best_qic),
                as.numeric(QIC(geem_afterph_wgt_list[[i]])),
                as.numeric(geem2_afterph_wgt_list[[i]]$best_qic),
                as.numeric(QIC(geem_afterph_alt_wgt_list[[i]])),
                as.numeric(geem2_afterph_alt_wgt_list[[i]]$best_qic),
                
                as.numeric(QIC(geem_strata_nowgt_list[[i]])),
                as.numeric(geem2_strata_nowgt_list[[i]]$best_qic),
                as.numeric(QIC(geem_strata_wgt_list[[i]])),
                as.numeric(geem2_strata_wgt_list[[i]]$best_qic),
                as.numeric(QIC(geem_strata_alt_wgt_list[[i]])),
                as.numeric(geem2_strata_alt_wgt_list[[i]]$best_qic),
                as.numeric(QIC(geem_strata_alt_alt_wgt_list[[i]])),
                as.numeric(geem2_strata_alt_alt_wgt_list[[i]]$best_qic)                        
         ),
         rbind.data.frame(
           dplyr::filter(summary2.geem(geem_nocorr_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem2_nocorr_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem_nocorr_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem2_nocorr_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem_nocorr_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem2_nocorr_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem_afterph_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem2_afterph_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem_afterph_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem2_afterph_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem_afterph_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem2_afterph_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
           
           dplyr::filter(summary2.geem(geem_strata_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem2_strata_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),                   
           dplyr::filter(summary2.geem(geem_strata_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem2_strata_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem_strata_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem2_strata_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem_strata_alt_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
           dplyr::filter(summary2.geem(geem2_strata_alt_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2")
         ))%>% 
  {
    copiar_nombres2(.)
    assign(paste0("result_data_geem_", gsub(" ","_",plan_names[i])), ., envir = .GlobalEnv)
    print(.)
  }


i=3
plan_names[i]
cbind.data.frame(model= types_of_model,
rbind.data.frame(
  dplyr::filter(tidy_geese_model(nocorr_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(nocorr_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(nocorr_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_alt_alt_wgt_list[[i]]), Term=="policonsumo2"))) %>% 
  {
    copiar_nombres2(.)  
    print(.)
  }
cbind.data.frame(model= types_of_model2,
           QIC= c(as.numeric(QIC(geem_nocorr_nowgt_list[[i]])),
                  as.numeric(geem2_nocorr_nowgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_nocorr_wgt_list[[i]])),
                  as.numeric(geem2_nocorr_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_nocorr_alt_wgt_list[[i]])),
                  as.numeric(geem2_nocorr_alt_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_afterph_nowgt_list[[i]])),
                  as.numeric(geem2_afterph_nowgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_afterph_wgt_list[[i]])),
                  as.numeric(geem2_afterph_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_afterph_alt_wgt_list[[i]])),
                  as.numeric(geem2_afterph_alt_wgt_list[[i]]$best_qic),
                  
                  as.numeric(QIC(geem_strata_nowgt_list[[i]])),
                  as.numeric(geem2_strata_nowgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_strata_wgt_list[[i]])),
                  as.numeric(geem2_strata_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_strata_alt_wgt_list[[i]])),
                  as.numeric(geem2_strata_alt_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_strata_alt_alt_wgt_list[[i]])),
                  as.numeric(geem2_strata_alt_alt_wgt_list[[i]]$best_qic)                        
           ),
           rbind.data.frame(
             dplyr::filter(summary2.geem(geem_nocorr_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_nocorr_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_nocorr_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_nocorr_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_nocorr_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_nocorr_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_afterph_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_afterph_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_afterph_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_afterph_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_afterph_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_afterph_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             
             dplyr::filter(summary2.geem(geem_strata_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),                   
             dplyr::filter(summary2.geem(geem_strata_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_strata_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_strata_alt_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_alt_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2")
           ))%>% 
  {
    copiar_nombres2(.)
    assign(paste0("result_data_geem_", gsub(" ","_",plan_names[i])), ., envir = .GlobalEnv)
    print(.)
  }


i=4
plan_names[i]
cbind.data.frame(model= types_of_model,
rbind.data.frame(
  dplyr::filter(tidy_geese_model(nocorr_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(nocorr_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(nocorr_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_alt_alt_wgt_list[[i]]), Term=="policonsumo2"))) %>% 
  {
    copiar_nombres2(.)  
    print(.)
  }
cbind.data.frame(model= types_of_model2,
           QIC= c(as.numeric(QIC(geem_nocorr_nowgt_list[[i]])),
                  as.numeric(geem2_nocorr_nowgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_nocorr_wgt_list[[i]])),
                  as.numeric(geem2_nocorr_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_nocorr_alt_wgt_list[[i]])),
                  as.numeric(geem2_nocorr_alt_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_afterph_nowgt_list[[i]])),
                  as.numeric(geem2_afterph_nowgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_afterph_wgt_list[[i]])),
                  as.numeric(geem2_afterph_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_afterph_alt_wgt_list[[i]])),
                  as.numeric(geem2_afterph_alt_wgt_list[[i]]$best_qic),
                  
                  as.numeric(QIC(geem_strata_nowgt_list[[i]])),
                  as.numeric(geem2_strata_nowgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_strata_wgt_list[[i]])),
                  as.numeric(geem2_strata_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_strata_alt_wgt_list[[i]])),
                  as.numeric(geem2_strata_alt_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_strata_alt_alt_wgt_list[[i]])),
                  as.numeric(geem2_strata_alt_alt_wgt_list[[i]]$best_qic)                        
           ),
           rbind.data.frame(
             dplyr::filter(summary2.geem(geem_nocorr_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_nocorr_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_nocorr_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_nocorr_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_nocorr_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_nocorr_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_afterph_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_afterph_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_afterph_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_afterph_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_afterph_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_afterph_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             
             dplyr::filter(summary2.geem(geem_strata_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),                   
             dplyr::filter(summary2.geem(geem_strata_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_strata_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_strata_alt_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_alt_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2")
           ))%>% 
  {
    copiar_nombres2(.)
    assign(paste0("result_data_geem_", gsub(" ","_",plan_names[i])), ., envir = .GlobalEnv)
    print(.)
  }

i=5
plan_names[i]
cbind.data.frame(model= types_of_model,
rbind.data.frame(
  dplyr::filter(tidy_geese_model(nocorr_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(nocorr_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(nocorr_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(afterph_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_nowgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_alt_wgt_list[[i]]), Term=="policonsumo2"),
  dplyr::filter(tidy_geese_model(strata_alt_alt_wgt_list[[i]]), Term=="policonsumo2"))) %>% 
  {
        copiar_nombres2(.)  
        print(.)
  }
  
cbind.data.frame(model= types_of_model2,
           QIC= c(as.numeric(QIC(geem_nocorr_nowgt_list[[i]])),
                  as.numeric(geem2_nocorr_nowgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_nocorr_wgt_list[[i]])),
                  as.numeric(geem2_nocorr_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_nocorr_alt_wgt_list[[i]])),
                  as.numeric(geem2_nocorr_alt_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_afterph_nowgt_list[[i]])),
                  as.numeric(geem2_afterph_nowgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_afterph_wgt_list[[i]])),
                  as.numeric(geem2_afterph_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_afterph_alt_wgt_list[[i]])),
                  as.numeric(geem2_afterph_alt_wgt_list[[i]]$best_qic),
                  
                  as.numeric(QIC(geem_strata_nowgt_list[[i]])),
                  as.numeric(geem2_strata_nowgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_strata_wgt_list[[i]])),
                  as.numeric(geem2_strata_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_strata_alt_wgt_list[[i]])),
                  as.numeric(geem2_strata_alt_wgt_list[[i]]$best_qic),
                  as.numeric(QIC(geem_strata_alt_alt_wgt_list[[i]])),
                  as.numeric(geem2_strata_alt_alt_wgt_list[[i]]$best_qic)                        
           ),
           rbind.data.frame(
             dplyr::filter(summary2.geem(geem_nocorr_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_nocorr_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_nocorr_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_nocorr_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_nocorr_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_nocorr_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_afterph_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_afterph_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_afterph_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_afterph_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_afterph_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_afterph_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             
             dplyr::filter(summary2.geem(geem_strata_nowgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_nowgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),                   
             dplyr::filter(summary2.geem(geem_strata_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_strata_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem_strata_alt_alt_wgt_list[[i]],exponentiate=T, digits=2),Term=="policonsumo2"),
             dplyr::filter(summary2.geem(geem2_strata_alt_alt_wgt_list[[i]]$best_model,exponentiate=T, digits=2),Term=="policonsumo2")
           ))%>% 
  {
    copiar_nombres2(.)
    assign(paste0("result_data_geem_", gsub(" ","_",plan_names[i])), ., envir = .GlobalEnv)
    print(.)
  }

Time for this code chunk to run: 0 minutes

Code
rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
                 cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
                 cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
                 ) %>% 
  dplyr::group_by(setting) %>% 
  dplyr::mutate(min_qic_setting=ifelse(QIC==min(QIC),1,0)) %>% 
  dplyr::mutate(QIC= format(as.numeric(sprintf("%1.1f",QIC)), big.mark=",")) %>% 
  dplyr::ungroup() %>%
  dplyr::mutate(`RR (95% IC)`= paste0(Estimates, " (", `CI.lower`, ", ", CI.upper, ")")) %>% 
  dplyr::select(setting, model, QIC, `RR (95% IC)`, p.value) %>% 
  {
  #copiar_nombres2()
    write.table(., file = paste0(folder_path,"gee_240515.csv"), dec=",", sep="\t")
    knitr::kable(., size=10, format="html", caption="Relative risk of treatment non-completion status by reported polysubstance use") %>% 
    kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>% 
    kableExtra::scroll_box(width = "100%", height = "375px")  
  }
Relative risk of treatment non-completion status by reported polysubstance use
setting model QIC RR (95% IC) p.value
basic ambulatory No time-varying coefficients, no weight 19,453.6 1.03 (1.00, 1.05) 0.0689
basic ambulatory No time-varying coefficients, no weight, NB 19,838.6 1.03 (1.00, 1.05) 0.0698
basic ambulatory No time-varying coefficients, lag=0 19,465.1 1.02 (0.99, 1.05) 0.1648
basic ambulatory No time-varying coefficients, lag=0, NB 19,465.1 1.02 (0.99, 1.05) 0.1648
basic ambulatory No time-varying coefficients, lag=1 19,458.8 1.02 (1.00, 1.05) 0.0779
basic ambulatory No time-varying coefficients, lag=1, NB 19,458.8 1.02 (1.00, 1.05) 0.0779
basic ambulatory With time-varying coefficients, no weight 19,453.6 1.03 (1.00, 1.05) 0.0689
basic ambulatory With time-varying coefficients, no weight, NB 19,453.6 1.03 (1.00, 1.05) 0.0689
basic ambulatory With time-varying coefficients, lag=0 19,479.4 1.03 (0.99, 1.06) 0.1328
basic ambulatory With time-varying coefficients, lag=0, NB 19,479.4 1.03 (0.99, 1.06) 0.1328
basic ambulatory With time-varying coefficients, lag=1 19,519.4 1.03 (0.98, 1.07) 0.2275
basic ambulatory With time-varying coefficients, lag=1, NB 19,519.4 1.03 (0.98, 1.07) 0.2275
basic ambulatory Stratified by follow-up intervals, no weight 19,453.6 1.03 (1.00, 1.05) 0.0689
basic ambulatory Stratified by follow-up intervals, no weight, NB 19,453.6 1.03 (1.00, 1.05) 0.0689
basic ambulatory Stratified by follow-up intervals, lag=0 19,476.8 1.02 (0.99, 1.06) 0.2682
basic ambulatory Stratified by follow-up intervals, lag=0, NB 19,476.8 1.02 (0.99, 1.06) 0.2682
basic ambulatory Stratified by follow-up intervals, lag=1 19,467.2 1.04 (1.01, 1.07) 0.0202
basic ambulatory Stratified by follow-up intervals, lag=1, NB 19,467.2 1.04 (1.01, 1.07) 0.0202
basic ambulatory Stratified by follow-up intervals, 2nd, lag=1 19,505.4 1.03 (0.99, 1.07) 0.1745
basic ambulatory Stratified by follow-up intervals, 2nd, lag=1, NB 19,505.4 1.03 (0.99, 1.07) 0.1745
GP intensive ambulatory No time-varying coefficients, no weight 22,261.9 1.04 (1.01, 1.07) 0.0082
GP intensive ambulatory No time-varying coefficients, no weight, NB 22,686.8 1.04 (1.01, 1.07) 0.008
GP intensive ambulatory No time-varying coefficients, lag=0 22,275.9 1.04 (1.01, 1.08) 0.01
GP intensive ambulatory No time-varying coefficients, lag=0, NB 22,276.0 1.04 (1.01, 1.08) 0.01
GP intensive ambulatory No time-varying coefficients, lag=1 22,271.9 1.04 (1.01, 1.07) 0.0113
GP intensive ambulatory No time-varying coefficients, lag=1, NB 22,271.9 1.04 (1.01, 1.07) 0.0113
GP intensive ambulatory With time-varying coefficients, no weight 22,261.9 1.04 (1.01, 1.07) 0.0082
GP intensive ambulatory With time-varying coefficients, no weight, NB 22,261.9 1.04 (1.01, 1.07) 0.0082
GP intensive ambulatory With time-varying coefficients, lag=0 22,284.4 1.06 (1.02, 1.09) 0.0032
GP intensive ambulatory With time-varying coefficients, lag=0, NB 22,284.4 1.06 (1.02, 1.09) 0.0032
GP intensive ambulatory With time-varying coefficients, lag=1 22,324.6 1.05 (1.00, 1.10) 0.0443
GP intensive ambulatory With time-varying coefficients, lag=1, NB 22,324.6 1.05 (1.00, 1.10) 0.0443
GP intensive ambulatory Stratified by follow-up intervals, no weight 22,261.9 1.04 (1.01, 1.07) 0.0082
GP intensive ambulatory Stratified by follow-up intervals, no weight, NB 22,261.9 1.04 (1.01, 1.07) 0.0082
GP intensive ambulatory Stratified by follow-up intervals, lag=0 22,283.0 1.04 (1.01, 1.08) 0.0174
GP intensive ambulatory Stratified by follow-up intervals, lag=0, NB 22,283.0 1.04 (1.01, 1.08) 0.0174
GP intensive ambulatory Stratified by follow-up intervals, lag=1 22,269.4 1.02 (0.99, 1.05) 0.1429
GP intensive ambulatory Stratified by follow-up intervals, lag=1, NB 22,269.4 1.02 (0.99, 1.05) 0.1429
GP intensive ambulatory Stratified by follow-up intervals, 2nd, lag=1 22,280.0 1.01 (0.98, 1.05) 0.4077
GP intensive ambulatory Stratified by follow-up intervals, 2nd, lag=1, NB 22,280.0 1.01 (0.98, 1.05) 0.4077
GP residential No time-varying coefficients, no weight 9,818.1 0.97 (0.92, 1.02) 0.2207
GP residential No time-varying coefficients, no weight, NB 9,979.2 0.97 (0.92, 1.02) 0.2203
GP residential No time-varying coefficients, lag=0 9,833.7 0.97 (0.92, 1.02) 0.2773
GP residential No time-varying coefficients, lag=0, NB 9,833.7 0.97 (0.92, 1.02) 0.2773
GP residential No time-varying coefficients, lag=1 9,827.8 0.95 (0.90, 1.01) 0.1236
GP residential No time-varying coefficients, lag=1, NB 9,827.8 0.95 (0.90, 1.01) 0.1236
GP residential With time-varying coefficients, no weight 9,818.1 0.97 (0.92, 1.02) 0.2207
GP residential With time-varying coefficients, no weight, NB 9,818.1 0.97 (0.92, 1.02) 0.2207
GP residential With time-varying coefficients, lag=0 9,830.3 0.98 (0.93, 1.04) 0.4692
GP residential With time-varying coefficients, lag=0, NB 9,830.3 0.98 (0.93, 1.04) 0.4692
GP residential With time-varying coefficients, lag=1 9,868.4 1.00 (0.94, 1.07) 0.991
GP residential With time-varying coefficients, lag=1, NB 9,868.4 1.00 (0.94, 1.07) 0.991
GP residential Stratified by follow-up intervals, no weight 9,818.1 0.97 (0.92, 1.02) 0.2207
GP residential Stratified by follow-up intervals, no weight, NB 9,818.1 0.97 (0.92, 1.02) 0.2207
GP residential Stratified by follow-up intervals, lag=0 9,839.4 0.95 (0.89, 1.02) 0.137
GP residential Stratified by follow-up intervals, lag=0, NB 9,839.4 0.95 (0.89, 1.02) 0.137
GP residential Stratified by follow-up intervals, lag=1 9,837.5 0.98 (0.91, 1.04) 0.4662
GP residential Stratified by follow-up intervals, lag=1, NB 9,837.5 0.98 (0.91, 1.04) 0.4662
GP residential Stratified by follow-up intervals, 2nd, lag=1 9,885.3 1.01 (0.92, 1.11) 0.8284
GP residential Stratified by follow-up intervals, 2nd, lag=1, NB 9,885.3 1.01 (0.92, 1.11) 0.8284
WO intensive ambulatory No time-varying coefficients, no weight 3,424.2 0.99 (0.93, 1.05) 0.715
WO intensive ambulatory No time-varying coefficients, no weight, NB 3,487.8 0.99 (0.92, 1.06) 0.7171
WO intensive ambulatory No time-varying coefficients, lag=0 3,434.7 0.99 (0.92, 1.07) 0.8848
WO intensive ambulatory No time-varying coefficients, lag=0, NB 3,434.7 0.99 (0.92, 1.07) 0.8848
WO intensive ambulatory No time-varying coefficients, lag=1 3,431.3 0.99 (0.92, 1.06) 0.7415
WO intensive ambulatory No time-varying coefficients, lag=1, NB 3,431.3 0.99 (0.92, 1.06) 0.7415
WO intensive ambulatory With time-varying coefficients, no weight 3,424.2 0.99 (0.93, 1.05) 0.715
WO intensive ambulatory With time-varying coefficients, no weight, NB 3,424.2 0.99 (0.93, 1.05) 0.715
WO intensive ambulatory With time-varying coefficients, lag=0 3,445.9 1.01 (0.92, 1.11) 0.8115
WO intensive ambulatory With time-varying coefficients, lag=0, NB 3,445.9 1.01 (0.92, 1.11) 0.8115
WO intensive ambulatory With time-varying coefficients, lag=1 3,468.6 1.01 (0.91, 1.12) 0.8504
WO intensive ambulatory With time-varying coefficients, lag=1, NB 3,468.6 1.01 (0.91, 1.12) 0.8504
WO intensive ambulatory Stratified by follow-up intervals, no weight 3,424.2 0.99 (0.93, 1.05) 0.715
WO intensive ambulatory Stratified by follow-up intervals, no weight, NB 3,424.2 0.99 (0.93, 1.05) 0.715
WO intensive ambulatory Stratified by follow-up intervals, lag=0 3,453.0 0.98 (0.90, 1.07) 0.6456
WO intensive ambulatory Stratified by follow-up intervals, lag=0, NB 3,453.0 0.98 (0.90, 1.07) 0.6456
WO intensive ambulatory Stratified by follow-up intervals, lag=1 3,437.3 0.98 (0.91, 1.05) 0.547
WO intensive ambulatory Stratified by follow-up intervals, lag=1, NB 3,437.3 0.98 (0.91, 1.05) 0.547
WO intensive ambulatory Stratified by follow-up intervals, 2nd, lag=1 3,466.2 0.94 (0.87, 1.03) 0.1744
WO intensive ambulatory Stratified by follow-up intervals, 2nd, lag=1, NB 3,466.2 0.94 (0.87, 1.03) 0.1744
WO residential No time-varying coefficients, no weight 4,826.1 1.14 (1.06, 1.23) 6e-04
WO residential No time-varying coefficients, no weight, NB 4,909.1 1.14 (1.06, 1.23) 6e-04
WO residential No time-varying coefficients, lag=0 4,839.1 1.15 (1.06, 1.25) 0.001
WO residential No time-varying coefficients, lag=0, NB 4,839.1 1.15 (1.06, 1.25) 0.001
WO residential No time-varying coefficients, lag=1 4,835.6 1.13 (1.04, 1.22) 0.0027
WO residential No time-varying coefficients, lag=1, NB 4,835.6 1.13 (1.04, 1.22) 0.0027
WO residential With time-varying coefficients, no weight 4,826.1 1.14 (1.06, 1.23) 6e-04
WO residential With time-varying coefficients, no weight, NB 4,826.1 1.14 (1.06, 1.23) 6e-04
WO residential With time-varying coefficients, lag=0 4,839.7 1.11 (1.02, 1.21) 0.0134
WO residential With time-varying coefficients, lag=0, NB 4,839.7 1.11 (1.02, 1.21) 0.0134
WO residential With time-varying coefficients, lag=1 4,861.0 1.13 (1.03, 1.25) 0.012
WO residential With time-varying coefficients, lag=1, NB 4,861.0 1.13 (1.03, 1.25) 0.012
WO residential Stratified by follow-up intervals, no weight 4,826.1 1.14 (1.06, 1.23) 6e-04
WO residential Stratified by follow-up intervals, no weight, NB 4,826.1 1.14 (1.06, 1.23) 6e-04
WO residential Stratified by follow-up intervals, lag=0 4,851.1 1.12 (1.03, 1.22) 0.0105
WO residential Stratified by follow-up intervals, lag=0, NB 4,851.1 1.12 (1.03, 1.22) 0.0105
WO residential Stratified by follow-up intervals, lag=1 4,841.8 1.11 (1.02, 1.20) 0.0117
WO residential Stratified by follow-up intervals, lag=1, NB 4,841.8 1.11 (1.02, 1.20) 0.0117
WO residential Stratified by follow-up intervals, 2nd, lag=1 4,874.2 1.09 (0.99, 1.20) 0.0693
WO residential Stratified by follow-up intervals, 2nd, lag=1, NB 4,874.2 1.09 (0.99, 1.20) 0.0693

Time for this code chunk to run: 0 minutes

Stratifying by presence of alcohol among poly-substance use

Code
data_mine_miss_restr_proc2exp2<-
  data_mine_miss_restr_proc2 %>% 
  left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","n_off_acq","n_off_vio","n_off_sud","n_off_oth", "n_prev_off")], by="hash_key") %>%
  tidylog::left_join(Base_fiscalia_v15f_grant_23_24[,c("hash_key","fech_ing_num","otras_sus1_mod","otras_sus2_mod","otras_sus3_mod")], by=c("hash_key"="hash_key","fech_ing_num"="fech_ing_num")) %>% 
  {if (nrow(.) > nrow(data_mine_miss_restr_proc2)) {message("left join added rows"); .} else .} %>% 
  dplyr::mutate(n_prev_off_bin= ifelse(n_prev_off>0,1,0)) %>% 
  dplyr::mutate(
    susinidumrec_coc = ifelse(sus_ini_mod_mvv == "Cocaine hydrochloride", 1, 0),
    susinidumrec_pbc = ifelse(sus_ini_mod_mvv == "Cocaine paste", 1, 0),
    susinidumrec_mar = ifelse(sus_ini_mod_mvv == "Marijuana", 1, 0),
    susinidumrec_otr = ifelse(sus_ini_mod_mvv == "Other", 1, 0), #sus_principal_mod
    susprindumrec_coc = ifelse(sus_principal_mod == "Cocaine hydrochloride", 1, 0),
    susprindumrec_pbc = ifelse(sus_principal_mod == "Cocaine paste", 1, 0),
    susprindumrec_mar = ifelse(sus_principal_mod == "Marijuana", 1, 0),
    susprindumrec_otr = ifelse(sus_principal_mod == "Other", 1, 0) # "susprindum_coc", "susprindum_pbc", "susprindum_mar", "susprindum_oh
  ) %>% 
  rowwise() %>%
  #Condition
  dplyr::mutate(oh_otras_sus = if_else("Alcohol" %in% c(otras_sus1_mod, otras_sus2_mod, otras_sus3_mod),1,0)) %>%
  dplyr::mutate(policonsumo2_rec= dplyr::case_when(policonsumo2==1 & oh_otras_sus==1~"2.both",policonsumo2==1 & oh_otras_sus==0~"1.only PSU",T~"0.no PSU")) %>% 
  ungroup()

left_join: added 3 columns (otras_sus1_mod, otras_sus2_mod, otras_sus3_mod)

       > rows only in x   18,026
       > rows only in y  (72,086)
       > matched rows     12,962
       >                 ========
       > rows total       30,988
Code
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#

nocorr_nowgt_int_list<-list()
nocorr_nowgt_int_probs_list<-list()
nocor_nowgt_int_tab <-list()
for (i in seq_along(plan_names)) {
  # Subset the data for the current plan type
  current_data <- subset(data_mine_miss_restr_proc2exp2, tipo_de_plan_2_mod == plan_names[i])
  
  # Fit the GEE model for the current subset
  model <- geeglm(tr_outcome ~ policonsumo2_rec +
                   comp_bpsc_y3_severe+
                   edad_al_ing_1 + 
                   ano_nac_corr + 
                   susinidumrec_otr +
                   susinidumrec_coc +
                   susinidumrec_pbc +
                   susinidumrec_mar +
                   psycom_dum_study +
                   psycom_dum_with +
                   freq_cons_dum_5day +
                   cond_oc_dum_2inact +
                   cond_oc_dum_3unemp +
                   susprindumrec_coc +
                   susprindumrec_pbc +
                   susprindumrec_mar +
                   susprindumrec_otr, #origen_ingreso_mod fis_comorbidity_icd_10 dg_trs_cons_sus_or
                 id = id, 
                 data = current_data,
                 family = poisson, 
                 corstr = "independence")

  # Assign the model to the list with a name based on the plan name
  model_name <- gsub(" ", "_", plan_names[i])
  model_name <- gsub("[^[:alnum:]_]", "", model_name)  # Clean up non-alphanumeric characters
  
  nocor_nowgt_int_tab[[paste0("tab_", model_name)]] <- table(current_data$policonsumo2_rec, current_data$tr_outcome) %>% data.frame()
  nocorr_nowgt_int_list[[paste("model", model_name, sep = "_")]] <- model
  emmeans_response <- emmeans(model, ~policonsumo2_rec,  rg.limit=2e5, type = "response")
  nocorr_nowgt_int_probs_list[[paste0("emmeans_response_int_", model_name)]] <- emmeans_response
  assign(paste0("emmeans_response_int_", model_name), emmeans_response, envir = .GlobalEnv)
}

nocor_nowgt_int_tab[1]
$tab_basic_ambulatory
        Var1 Var2 Freq
1   0.no PSU    0  629
2 1.only PSU    0 1210
3     2.both    0  286
4   0.no PSU    1 1827
5 1.only PSU    1 4450
6     2.both    1 1591
Code
nocor_nowgt_int_tab[2]
$tab_GP_intensive_ambulatory
        Var1 Var2 Freq
1   0.no PSU    0  624
2 1.only PSU    0 1615
3     2.both    0  401
4   0.no PSU    1 1573
5 1.only PSU    1 5226
6     2.both    1 2058
Code
nocor_nowgt_int_tab[3]
$tab_GP_residential
        Var1 Var2 Freq
1   0.no PSU    0  261
2 1.only PSU    0  905
3     2.both    0  400
4   0.no PSU    1  572
5 1.only PSU    1 2360
6     2.both    1  706
Code
nocor_nowgt_int_tab[4]
$tab_WO_intensive_ambulatory
        Var1 Var2 Freq
1   0.no PSU    0  110
2 1.only PSU    0  253
3     2.both    0   61
4   0.no PSU    1  292
5 1.only PSU    1  726
6     2.both    1  318
Code
nocor_nowgt_int_tab[5]
$tab_WO_residential
        Var1 Var2 Freq
1   0.no PSU    0  172
2 1.only PSU    0  413
3     2.both    0  130
4   0.no PSU    1  289
5 1.only PSU    1 1147
6     2.both    1  383
Code
# > nocor_nowgt_int_tab[1]
# $tab_basic_ambulatory
# Var1 Var2 Freq
# 1   0.no PSU    0  629
# 2 1.only PSU    0 1210
# 3     2.both    0  286
# 4   0.no PSU    1 1827
# 5 1.only PSU    1 4450
# 6     2.both    1 1591
# 
# > nocor_nowgt_int_tab[2]
# $tab_GP_intensive_ambulatory
# Var1 Var2 Freq
# 1   0.no PSU    0  624
# 2 1.only PSU    0 1615
# 3     2.both    0  401
# 4   0.no PSU    1 1573
# 5 1.only PSU    1 5226
# 6     2.both    1 2058
# 
# > nocor_nowgt_int_tab[3]
# $tab_GP_residential
# Var1 Var2 Freq
# 1   0.no PSU    0  261
# 2 1.only PSU    0  905
# 3     2.both    0  400
# 4   0.no PSU    1  572
# 5 1.only PSU    1 2360
# 6     2.both    1  706
# 
# > nocor_nowgt_int_tab[4]
# $tab_WO_intensive_ambulatory
# Var1 Var2 Freq
# 1   0.no PSU    0  110
# 2 1.only PSU    0  253
# 3     2.both    0   61
# 4   0.no PSU    1  292
# 5 1.only PSU    1  726
# 6     2.both    1  318
# 
# > nocor_nowgt_int_tab[5]
# $tab_WO_residential
# Var1 Var2 Freq
# 1   0.no PSU    0  172
# 2 1.only PSU    0  413
# 3     2.both    0  130
# 4   0.no PSU    1  289
# 5 1.only PSU    1 1147
# 6     2.both    1  383
# 
invisible("no hay tan pocos casos en un estrato, salvo en WO intensive ambulatory")



invisible("En el basic ambulatory, GP intensive-amb y GP residential, el sig es el ambos (consumo alcohol)")
  
#broom::tidy(nocorr_nowgt_int_list[[1]], exponentiate=T, conf.int=T)[1:3,] 
#    term                         estimate std.error statistic     p.value conf.low conf.high
#   <chr>                           <dbl>     <dbl>     <dbl>       <dbl>    <dbl>     <dbl>
# 1 (Intercept)                0.00000870    4.83       5.81  0.0159      6.68e-10     0.113
# 2 policonsumo2_rec1.only PSU 1.01          0.0143     0.488 0.485       9.82e- 1     1.04 
# 3 policonsumo2_rec2.both     1.08          0.0161    24.8   0.000000622 1.05e+ 0     1.12 

#broom::tidy(nocorr_nowgt_int_list[[2]], exponentiate=T, conf.int=T)[1:3,]
#   term                         estimate std.error statistic       p.value conf.low conf.high
#   <chr>                           <dbl>     <dbl>     <dbl>         <dbl>    <dbl>     <dbl>
# 1 (Intercept)                0.00000655    4.81        6.15 0.0132        5.23e-10    0.0822
# 2 policonsumo2_rec1.only PSU 1.02          0.0153      1.98 0.160         9.92e- 1    1.05  
# 3 policonsumo2_rec2.both     1.10          0.0162     35.1  0.00000000312 1.07e+ 0    1.14  

invisible("Aunque GP residential es protector el policonsumo en both")
#broom::tidy(nocorr_nowgt_int_list[[3]], exponentiate=T, conf.int=T)[1:3,]
#   term                       estimate std.error statistic  p.value     conf.low    conf.high
#   <chr>                         <dbl>     <dbl>     <dbl>    <dbl>        <dbl>        <dbl>
# 1 (Intercept)                   1.14     9.06    0.000199 0.989    0.0000000221 58484383.   
# 2 policonsumo2_rec1.only PSU    0.989    0.0275  0.175    0.676    0.937               1.04 
# 3 policonsumo2_rec2.both        0.885    0.0332 13.5      0.000243 0.830               0.945


invisible("En WO intensive ambulatory, no tira ninguno")

#broom::tidy(nocorr_nowgt_int_list[[4]], exponentiate=T, conf.int=T)[1:3,]
# # A tibble: 3 x 7
#   term                       estimate std.error statistic p.value conf.low  conf.high
#   <chr>                         <dbl>     <dbl>     <dbl>   <dbl>    <dbl>      <dbl>
# 1 (Intercept)                0.000318   12.0        0.452  0.501  2.02e-14 4991312.  
# 2 policonsumo2_rec1.only PSU 0.956       0.0357     1.62   0.203  8.91e- 1       1.02
# 3 policonsumo2_rec2.both     1.07        0.0376     3.13   0.0767 9.93e- 1       1.15
# 
invisible("En WO residential, ambos mueven la aguja a mayor riesgo")
# 
#broom::tidy(nocorr_nowgt_int_list[[5]], exponentiate=T, conf.int=T)[1:3,]
#   term                        estimate std.error statistic  p.value  conf.low conf.high
#   <chr>                          <dbl>     <dbl>     <dbl>    <dbl>     <dbl>     <dbl>
# 1 (Intercept)                249351.     11.6         1.15 0.284    0.0000331   1.88e15
# 2 policonsumo2_rec1.only PSU      1.14    0.0388     11.3  0.000775 1.06        1.23e 0
# 3 policonsumo2_rec2.both          1.14    0.0439      8.83 0.00297  1.05        1.24e 0
# 
# 
bind_rows(
    lapply(nocorr_nowgt_int_list[1:5], function(model) {
        broom::tidy(model, exponentiate = TRUE, conf.int = TRUE) %>%
            dplyr::select(term, estimate, conf.low, conf.high, p.value) %>%
        
            slice(2:3)
    }),
    .id = "model"
) %>% 
  dplyr::mutate(term= gsub("policonsumo2_","", term)) %>% 
  dplyr::mutate(across(c("estimate","conf.low", "conf.high"),~sprintf("%1.2f",.))) %>% 
  dplyr::mutate(RR= paste0(estimate, " (", conf.low,"-", conf.high,")")) %>% 
  dplyr::select(-estimate, -conf.low, -conf.high) %>% 
  dplyr::select(model, term, RR, p.value) %>% 
  dplyr::mutate(across(c("p.value"),~sprintf("%1.4f",.))) %>%
  {
  knitr::kable(., "markdown", caption= "Relative risk of treatment non-completion status by reported polysubstance use, whether ")
  }
Relative risk of treatment non-completion status by reported polysubstance use, whether
model term RR p.value
model_basic_ambulatory rec1.only PSU 1.01 (0.98-1.04) 0.4849
model_basic_ambulatory rec2.both 1.08 (1.05-1.12) 0.0000
model_GP_intensive_ambulatory rec1.only PSU 1.02 (0.99-1.05) 0.1599
model_GP_intensive_ambulatory rec2.both 1.10 (1.07-1.14) 0.0000
model_GP_residential rec1.only PSU 0.99 (0.94-1.04) 0.6761
model_GP_residential rec2.both 0.89 (0.83-0.94) 0.0002
model_WO_intensive_ambulatory rec1.only PSU 0.96 (0.89-1.02) 0.2034
model_WO_intensive_ambulatory rec2.both 1.07 (0.99-1.15) 0.0767
model_WO_residential rec1.only PSU 1.14 (1.06-1.23) 0.0008
model_WO_residential rec2.both 1.14 (1.05-1.24) 0.0030

Time for this code chunk to run: 0.3 minutes

Code
# bind_rows(
#   lapply(1:5, function(model) {
#     nocor_nowgt_int_tab[[model]]%>% dplyr::mutate(Var0=ifelse(grepl("no PSU",Var1),1,0)) %>% dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n= sum(Freq))
#   }),
#   .id = "model"
# ) %>% 
#   group_by(model) %>% 
#   summarise(n=sum(n))
data_barplot <- bind_rows(
  lapply(1:5, function(model) {
    nocor_nowgt_int_tab[[model]] %>% 
      dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% 
      dplyr::group_by(Var0, Var2) %>% 
      dplyr::summarise(n = sum(Freq))
  }),
  .id = "model"
) %>% 
  group_by(model, Var0) %>%
  mutate(percentage = n / sum(n)) %>%
  ungroup()

summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument.

Code
# Filter data for treatment non-completion only
data_non_completion <- data_barplot %>% filter(Var2 == "1")

# Create labels for the total counts and percentages for each group
labels <- data_barplot %>%
  group_by(model, Var0) %>%
  summarise(total = sum(n)) %>%
  mutate(label = scales::percent(total / sum(total), accuracy = 0.1))

summarise() has grouped output by ‘model’. You can override using the .groups argument.

Code
# Calculate overall percentages for PSU and no PSU
overall_percentages <- data_barplot %>%
  group_by(model, Var0) %>%
  summarise(total = sum(n)) %>%
  mutate(overall_percentage = total / sum(total))

summarise() has grouped output by ‘model’. You can override using the .groups argument.

Code
# Create a separate data frame for the labels
labels <- overall_percentages %>%
  mutate(label = scales::percent(overall_percentage, accuracy = 1))

data_barplot$Var0 <- factor(data_barplot$Var0, levels = c("Single\nsubstance use", "Polysubstance\nuse"))
labels$Var0 <- factor(labels$Var0, levels = c("Single\nsubstance use", "Polysubstance\nuse"))

showtext_auto()
showtext_opts(dpi = 500)
invisible("Si no funciona el times new roman")
#windowsFonts(TimesNewRoman = windowsFont("Times New Roman"))
# names(grDevices::windowsFonts())
# Generate the plot
ggplot(data_barplot, aes(x = factor(Var0), y = n, fill = factor(Var2))) +
  geom_bar(stat = "identity", position = position_stack(reverse = TRUE)) +
  #geom_bar(stat = "identity", position = "stack") +
  facet_wrap(~ model, ncol = 1, labeller = labeller(model = c(
    "1" = "Basic ambulatory (n= 9,993)",
    "2" = "General-population, intensive ambulatory (n= 11,497)",
    "3" = "General-population, residential (n= 5,204)",
    "4" = "Women-specific, intensive ambulatory (n= 1,760)",
    "5" = "Women-specific, residential (n= 2,534)"
  )), strip.position = "top") +
  scale_fill_manual(values = c("grey80", "grey35"), label=c("Completion", "Non-completion")) +
  labs(x = "Polysubstance use (PSU)", y = "Count", fill = "Treatment\noutcome") +
  geom_label(data = data_non_completion, aes(label = scales::percent(percentage, accuracy = 1)), position = position_stack(vjust = 0.8, reverse = TRUE), size = 3.9, fontface = "bold", color = "white", fill = "black", label.size = 0.2, family = "serif") +
   geom_text(data = labels, aes(x = Var0, y = 9e3, label = label, fill="1"), 
             size = 3.9, color = "black", hjust = 0.5, vjust = -0.5, fontface = "bold.italic", family = "serif") +
  ylim(0,1.1e4)+
  #theme_minimal(base_size = 14, base_family = "Times New Roman") +
  theme(
    strip.text = element_text(size = 12, face = "bold", family = "serif"),
    #axis.text.x = element_text(angle = 0, hjust = 1),
    axis.title.x = element_blank(),
    axis.text = element_text(size = 12, family = "serif"),
    axis.text.x = element_text(family = "serif"),
    text = element_text(family = "serif"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    strip.text.x = element_text(angle = 0, hjust = 0.5, size = 12, face = "bold", family = "serif"),
    strip.background = element_blank(),
    axis.line = element_line(colour = "black"),
    legend.position = "bottom",
    legend.text = element_text(size= 12, family = "serif"),
    legend.title = element_text(size= 13, family = "serif"),
  )

Warning in geom_text(data = labels, aes(x = Var0, y = 9000, label = label, : Ignoring unknown aesthetics: fill

Code
ggsave("E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/_figs/barplot3.png", height=14*.68, width=10*.68, dpi=500)
ggsave("E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/_figs/barplot3.jpg", height=14*.68, width=10*.68, dpi=500)
Figure 1: Percentage distribution of treatment outcomes by treatment setting and polysubstance use status

Time for this code chunk to run: 0.2 minutes

Code
#data_barplot
#
mantelhaen.test(
xtabs(n ~Var2+ Var0 +model, data= data_barplot)
)

Warning in (x[1L, 1L, ] + x[2L, 2L, ]) * x[1L, 1L, ] * x[2L, 2L, ]: NAs producidos por enteros excedidos

Warning in (x[1L, 1L, ] + x[2L, 2L, ]) * x[1L, 2L, ] * x[2L, 1L, ]: NAs producidos por enteros excedidos

Warning in (x[1L, 2L, ] + x[2L, 1L, ]) * x[1L, 1L, ] * x[2L, 2L, ]: NAs producidos por enteros excedidos

Warning in (x[1L, 2L, ] + x[2L, 1L, ]) * x[1L, 2L, ] * x[2L, 1L, ]: NAs producidos por enteros excedidos


    Mantel-Haenszel chi-squared test with continuity correction

data:  xtabs(n ~ Var2 + Var0 + model, data = data_barplot)
Mantel-Haenszel X-squared = 93.675, df = 1, p-value < 2.2e-16
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 NA NA
sample estimates:
common odds ratio 
         1.361499 
Code
#Mantel-Haenszel X-squared = 93.675, df = 1, p-value < 2.2e-16

invisible("Si bien se ve una aso")

chisq.test(xtabs(n ~Var2+ Var0, data= data_barplot))

    Pearson's Chi-squared test with Yates' continuity correction

data:  xtabs(n ~ Var2 + Var0, data = data_barplot)
X-squared = 76.039, df = 1, p-value < 2.2e-16
Code
#X-squared = 76.039, df = 1, p-value < 2.2e-16

vcd::woolf_test(
  xtabs(n ~Var2+ Var0 +model, data= data_barplot)
)

    Woolf-test on Homogeneity of Odds Ratios (no 3-Way assoc.)

data:  xtabs(n ~ Var2 + Var0 + model, data = data_barplot)
X-squared = 13.735, df = 4, p-value = 0.008192
Code
#Woolf-test on Homogeneity of Odds Ratios (no 3-Way assoc.)
#X-squared = 13.735, df = 4, p-value = 0.008192

invisible("Los odds de no completar varían por baseline tr. setting (p<0,001)")
invisible("No se puede usar Maentel-haenzel")
invisible("Suggests conditional dependence and varying strength or direction of association across strata (interaction), cautioning against a simple summary of the association.[common OR may not be reliable]")


# nocor_nowgt_int_tab[[1]] %>% 
#       dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% 
#       dplyr::group_by(Var0, Var2) %>% 
#       dplyr::summarise(n = sum(Freq))
chisq.test(
xtabs(n ~ Var2+Var0, data = nocor_nowgt_int_tab[[1]] %>% 
          dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% 
          dplyr::group_by(Var0, Var2) %>% 
          dplyr::summarise(n = sum(Freq)))
)

summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument.


    Pearson's Chi-squared test with Yates' continuity correction

data:  xtabs(n ~ Var2 + Var0, data = nocor_nowgt_int_tab[[1]] %>% dplyr::mutate(Var0 = ifelse(grepl("no PSU",     Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%     dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n = sum(Freq)))
X-squared = 36.389, df = 1, p-value = 1.616e-09
Code
#X-squared = 36.389, df = 1, p-value = 1.616e-09
# Phi-Coefficient   : 0.061 
# Contingency Coeff.: 0.061 
# Cramer's V        : 0.061 


chisq.test(
xtabs(n ~ Var2+Var0, data = nocor_nowgt_int_tab[[2]] %>% 
          dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% 
          dplyr::group_by(Var0, Var2) %>% 
          dplyr::summarise(n = sum(Freq)))
)

summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument.


    Pearson's Chi-squared test with Yates' continuity correction

data:  xtabs(n ~ Var2 + Var0, data = nocor_nowgt_int_tab[[2]] %>% dplyr::mutate(Var0 = ifelse(grepl("no PSU",     Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%     dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n = sum(Freq)))
X-squared = 45.055, df = 1, p-value = 1.916e-11
Code
#X-squared = 45.055, df = 1, p-value = 1.916e-11
# Phi-Coefficient   : 0.063 
# Contingency Coeff.: 0.063 
# Cramer's V        : 0.063

chisq.test(
xtabs(n ~ Var2+Var0, data = nocor_nowgt_int_tab[[3]] %>% 
          dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% 
          dplyr::group_by(Var0, Var2) %>% 
          dplyr::summarise(n = sum(Freq)))
)

summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument.


    Pearson's Chi-squared test with Yates' continuity correction

data:  xtabs(n ~ Var2 + Var0, data = nocor_nowgt_int_tab[[3]] %>% dplyr::mutate(Var0 = ifelse(grepl("no PSU",     Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%     dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n = sum(Freq)))
X-squared = 0.65673, df = 1, p-value = 0.4177
Code
#X-squared = 0.65673, df = 1, p-value = 0.4177
# hi-Coefficient   : 0.012 
# Contingency Coeff.: 0.012 
# Cramer's V        : 0.012 
# 
chisq.test(
xtabs(n ~ Var2+Var0, data = nocor_nowgt_int_tab[[4]] %>% 
          dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% 
          dplyr::group_by(Var0, Var2) %>% 
          dplyr::summarise(n = sum(Freq)))
)

summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument.


    Pearson's Chi-squared test with Yates' continuity correction

data:  xtabs(n ~ Var2 + Var0, data = nocor_nowgt_int_tab[[4]] %>% dplyr::mutate(Var0 = ifelse(grepl("no PSU",     Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%     dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n = sum(Freq)))
X-squared = 2.8231, df = 1, p-value = 0.09291
Code
#X-squared = 2.8231, df = 1, p-value = 0.09291
# Phi-Coefficient   : 0.042 
# Contingency Coeff.: 0.042 
# Cramer's V        : 0.042 
#
invisible("setting 5= No completar vs. completar, según PSU vs. no-PSU")
chisq.test(
xtabs(n ~ Var2+Var0, data = nocor_nowgt_int_tab[[5]] %>% 
          dplyr::mutate(Var0 = ifelse(grepl("no PSU", Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>% 
          dplyr::group_by(Var0, Var2) %>% 
          dplyr::summarise(n = sum(Freq)))
)

summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument. summarise() has grouped output by ‘Var0’. You can override using the .groups argument.


    Pearson's Chi-squared test with Yates' continuity correction

data:  xtabs(n ~ Var2 + Var0, data = nocor_nowgt_int_tab[[5]] %>% dplyr::mutate(Var0 = ifelse(grepl("no PSU",     Var1), "Single\nsubstance use", "Polysubstance\nuse")) %>%     dplyr::group_by(Var0, Var2) %>% dplyr::summarise(n = sum(Freq)))
X-squared = 22.463, df = 1, p-value = 2.142e-06
Code
#X-squared = 22.463, df = 1, p-value = 2.142e-06
# Phi-Coefficient   : 0.095 
# Contingency Coeff.: 0.095 
# Cramer's V        : 0.095 
invisible("compare setting vs. PSUs")
chisq.test(xtabs(n ~Var0+ model, data= data_barplot))

    Pearson's Chi-squared test

data:  xtabs(n ~ Var0 + model, data = data_barplot)
X-squared = 194.31, df = 4, p-value < 2.2e-16
Code
#X-squared = 194.31, df = 4, p-value < 2.2e-16


# Cox, M. K., & Key, C. H. (1993). Post hoc pair-wise comparisons for the chi-square test of homogeneity of proportions. Educational and Psychological Measurement, 53(4), 951–962. https://doi.org/10.1177/0013164493053004008

invisible("Settings vs. noncompletion")

chisq.test(xtabs(n ~Var2+ model, data= data_barplot))

    Pearson's Chi-squared test

data:  xtabs(n ~ Var2 + model, data = data_barplot)
X-squared = 177.64, df = 4, p-value < 2.2e-16
Code
#X-squared = 177.64, df = 4, p-value < 2.2e-16

cat("For each baseline setting, calculate the std. residuals of chi-square, to compare non-completion shares")
For each baseline setting, calculate the std. residuals of chi-square, to compare non-completion shares
Code
data.frame(chisq.posthoc.test(xtabs(n ~Var2+ model, data= data_barplot))) |> 
    dplyr::mutate(X4=as.character(X4)) |> 
    pivot_longer(cols = c(`X1`,`X2`,`X3`,`X4`,`X5`), names_to = "Group", values_to = "Values") %>%
    pivot_wider(names_from = Value, values_from = Values) %>%
    dplyr::filter(Dimension=="1") %>%
  dplyr::mutate(across(1, ~ gsub("\n"," ",as.character(.))))|>
    dplyr::mutate(Residuals= sprintf("%1.2f", as.numeric(Residuals))) %>% 
    dplyr::mutate( `p values`= sprintf("%1.3f", as.numeric(gsub("\\*","",`p values`)))) %>%
    dplyr::select(Dimension, Residuals, `p values`) %>% 
  dplyr::mutate(`p values`= ifelse(`p values`=="0.000","<0.001",`p values`))
# A tibble: 5 x 3
  Dimension Residuals `p values`
  <chr>     <chr>     <chr>     
1 1         8.07      <0.001    
2 1         3.61      0.003     
3 1         -11.07    <0.001    
4 1         0.02      1.000     
5 1         -5.05     <0.001    
Code
data.frame(chisq.posthoc.test(xtabs(n ~Var0+ model, data= data_barplot))) |> 
    dplyr::mutate(X4=as.character(X4)) |> 
    pivot_longer(cols = c(`X1`,`X2`,`X3`,`X4`,`X5`), names_to = "Group", values_to = "Values") %>%
    pivot_wider(names_from = Value, values_from = Values) %>%
    dplyr::filter(Dimension=="Polysubstance\nuse") %>%
  dplyr::mutate(across(1, ~ gsub("\n"," ",as.character(.))))|>
    dplyr::mutate(Residuals= sprintf("%1.2f", as.numeric(Residuals))) %>% 
    dplyr::mutate( `p values`= sprintf("%1.3f", as.numeric(gsub("\\*","",`p values`)))) %>%
    dplyr::select(Dimension, Residuals, `p values`) %>% 
  dplyr::mutate(`p values`= ifelse(`p values`=="0.000","<0.001",`p values`))
# A tibble: 5 x 3
  Dimension         Residuals `p values`
  <chr>             <chr>     <chr>     
1 Polysubstance use -12.30    <0.001    
2 Polysubstance use 4.62      <0.001    
3 Polysubstance use 8.78      <0.001    
4 Polysubstance use -2.52     0.118     
5 Polysubstance use 2.99      0.028     

Time for this code chunk to run: 0 minutes

Code
#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#
nocorr_nowgt_int2_list<-
bind_rows(
  lapply(nocorr_nowgt_int_list[1:5], function(model) {
    broom::tidy(model, exponentiate = TRUE, conf.int = TRUE) %>%
      dplyr::select(term, estimate, std.error, conf.low, conf.high, p.value) %>%
      dplyr::mutate(log_rr= log(estimate), log_lcl= log(conf.low), log_ucl= log(conf.high), variance= ((log_ucl-log_rr)/qnorm(0.975))^2)%>% 
      dplyr::slice(2:3)
  }),
  .id = "model"
)

log_rr_int<- nocorr_nowgt_int2_list$log_rr  
variances_int<- nocorr_nowgt_int2_list$variance  


results_int <- expand.grid(i=1:5, j=1:5) %>%
  filter(i < j) %>%
  mutate(
    drrr = log_rr_int[i] - log_rr_int[j],
    var_drrr = variances_int[i] + variances_int[j],
    se_drrr = sqrt(var_drrr),
    z_value = drrr / se_drrr,
    p_value = 2 * (1 - pnorm(abs(z_value))),
    ratio_estimates = exp(drrr),
    lower_bound = exp(drrr - qnorm(0.975) * se_drrr),
    upper_bound = exp(drrr + qnorm(0.975) * se_drrr)
  ) %>% 
  dplyr::mutate(i= dplyr::case_when(i==1~plan_names[1], i==2~plan_names[2], 
                i==3~plan_names[3], i==4~plan_names[4], i==5~plan_names[5]), 
                j= dplyr::case_when(j==1~plan_names[1], j==2~plan_names[2], 
                j==3~plan_names[3], j==4~plan_names[4], j==5~plan_names[5]))

Time for this code chunk to run: 0 minutes

Marginal probabilities

Loop over each unique plan name to fit a model. We used an interaction of PSU and treatment setting to capture the heterogeneity of PSU by treatment setting.

Code
df_marginal_plots<-
  cbind.data.frame(wgt= c("wgt", "iiw_nocorr_st", "iiw_nocorr_alt_st", "wgt", "iiw_after_ph_st", "iiw_after_ph_st"), 
        df= c("data_mine_miss_restr_proc2exp","data_mine_miss_restr_proc2exp","data_mine_miss_restr_proc2exp","data_mine_miss_restr_proc2exp_iiw_after_ph","data_mine_miss_restr_proc2exp_iiw_after_ph","data_mine_miss_restr_proc2exp_iiw_after_ph_alt"),
        model= c("no corr no wgt","no corr main","no corr alt", "after ph no wgt", "after ph main", "after ph alt"))

# Loop over each unique plan name to fit a model
for (i in seq_along(df_marginal_plots$wgt)) {
  # Subset the data for the current plan type
  #current_data <- subset(data_mine_miss_restr_proc2, tipo_de_plan_2_mod == plan_names[i])
  geeglm( tr_outcome  ~ policonsumo2 * tipo_de_plan_2_mod +
                        edad_al_ing_1 + 
                        ano_nac_corr + 
                        susinidumrec_otr +
                        susinidumrec_coc +
                        susinidumrec_pbc +
                        susinidumrec_mar +
                        psycom_dum_study +
                        psycom_dum_with +
                        freq_cons_dum_5day +
                        cond_oc_dum_2inact +
                        cond_oc_dum_3unemp +
                       susprindumrec_coc +
                       susprindumrec_pbc +
                       susprindumrec_mar +
                       susprindumrec_otr, 
                    id= id, 
                    data= get(df_marginal_plots$df[i]) %>% dplyr::mutate(wgt=1), 
                    weight= get(df_marginal_plots$df[i]) %>% dplyr::mutate(wgt=1) %>% pull(df_marginal_plots$wgt[i]), 
                    family= poisson, 
                    corstr= "independence") %>% 
  assign(paste0("geeglm_", gsub(" ","_",df_marginal_plots$model[i])), ., envir = .GlobalEnv)
  print(message(paste0("Models: ",paste0("geeglm_", gsub(" ","_",df_marginal_plots$model[i])))))
  emmeans_response <- emmeans(get(paste0("geeglm_", gsub(" ","_",df_marginal_plots$model[i]))), ~policonsumo2+tipo_de_plan_2_mod,  rg.limit=2e5, type = "response")%>% assign(paste0("emmeans_response_", gsub(" ","_",df_marginal_plots$model[i])), ., envir = .GlobalEnv)
  print(message(paste0("Emmeans: ",paste0("emmeans_response_", gsub(" ","_",df_marginal_plots$model[i])))))
}
NULL
NULL
NULL
NULL
NULL
NULL
NULL
NULL
NULL
NULL
NULL
NULL

Time for this code chunk to run: 0.5 minutes

Code
plot(emmeans_response_no_corr_no_wgt)+ ggtitle("No corr, no weights")#+ geom_point(aes(x = percent, y = conc), data = pigs, pch = 2, color = "blue")

Code
plot(emmeans_response_no_corr_main)+ ggtitle("No corr, main weights")

Code
plot(emmeans_response_no_corr_alt)+ ggtitle("No corr, alt weights")

Code
plot(emmeans_response_after_ph_no_wgt)+ ggtitle("After PH, no weights")

Code
plot(emmeans_response_after_ph_main)+ ggtitle("After PH, main weights")

Code
plot(emmeans_response_after_ph_alt)+ ggtitle("After PH, alt weights")

Time for this code chunk to run: 0.6 minutes

Heterogeneity

Code
library(meta)
try(library(dmetar))
Error in library(dmetar) : there is no package called 'dmetar'
Code
library(metafor)


#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##
# Define los valores de los Relative Risk Reductions (RRR) y sus intervalos de confianza
rrr <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
                 cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
                 cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>% 
    dplyr::group_by(setting) %>% 
    dplyr::filter(model=="No time-varying coefficients, no weight") %>% pull(`Estimates`) %>% as.numeric()
lcl <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
                 cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
                 cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>% 
    dplyr::group_by(setting) %>% 
    dplyr::filter(model=="No time-varying coefficients, no weight") %>% pull(`CI.lower`) %>% as.numeric()
ucl <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
                 cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
                 cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>% 
    dplyr::group_by(setting) %>% 
    dplyr::filter(model=="No time-varying coefficients, no weight") %>% pull(`CI.upper`) %>% as.numeric()


# Calcula los logaritmos de los RRR y de los límites de sus intervalos de confianza
log_rrr <- log(rrr)
log_lcl <- log(lcl)
log_ucl <- log(ucl)

# Calcula las varianzas
variances <- ((log_ucl - log_rrr) / qnorm(0.975))^2

# Meta-análisis usando un modelo de efectos aleatorios
meta_analysis <- rma(yi=log_rrr, sei=sqrt(variances), method="REML")

# Prueba de Cochran Q para evaluar la heterogeneidad
cochrans_q <- meta_analysis$QE
p_value_q <- meta_analysis$QEp

cat(sprintf("No time-varying coefficients, no weight- Cochran's Q = %.2f, p-value = %.4f\n", cochrans_q, p_value_q))
No time-varying coefficients, no weight- Cochran's Q = 14.49, p-value = 0.0059
Code
#concuerda con el calculo de stata
#
#

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
###_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
###
###

# Define los valores de los Relative Risk Reductions (RRR) y sus intervalos de confianza
rrr2 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
                         cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
                         cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
                         cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
                         cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>% 
    dplyr::group_by(setting) %>% 
    dplyr::filter(model=="No time-varying coefficients, lag=0") %>% pull(`Estimates`) %>% as.numeric()
lcl2 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
                         cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
                         cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
                         cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
                         cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>% 
    dplyr::group_by(setting) %>% 
    dplyr::filter(model=="No time-varying coefficients, lag=0") %>% pull(`CI.lower`) %>% as.numeric()
ucl2 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
                         cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
                         cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
                         cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
                         cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>% 
    dplyr::group_by(setting) %>% 
    dplyr::filter(model=="No time-varying coefficients, lag=0") %>% pull(`CI.upper`) %>% as.numeric()


# Calcula los logaritmos de los RRR y de los límites de sus intervalos de confianza
log_rrr2 <- log(rrr2)
log_lcl2 <- log(lcl2)
log_ucl2 <- log(ucl2)

# Calcula las varianzas
variances2 <- ((log_ucl2 - log_rrr2) / qnorm(0.975))^2

# Meta-análisis usando un modelo de efectos aleatorios
meta_analysis2 <- rma(yi=log_rrr2, sei=sqrt(variances2), method="REML")

# Prueba de Cochran Q para evaluar la heterogeneidad
cochrans_q2 <- meta_analysis2$QE
p_value_q2 <- meta_analysis2$QEp

cat(sprintf("No time-varying coefficients, lag=0, main- Cochran's Q = %.2f, p-value = %.4f\n", cochrans_q2, p_value_q2))
No time-varying coefficients, lag=0, main- Cochran's Q = 13.32, p-value = 0.0098
Code
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
###_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

# Define los valores de los Relative Risk Reductions (RRR) y sus intervalos de confianza
rrr3 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
                 cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
                 cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>% 
    dplyr::group_by(setting) %>% 
    dplyr::filter(model=="No time-varying coefficients, lag=1") %>% pull(`Estimates`) %>% as.numeric()
lcl3 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
                 cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
                 cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>% 
    dplyr::group_by(setting) %>% 
    dplyr::filter(model=="No time-varying coefficients, lag=1") %>% pull(`CI.lower`) %>% as.numeric()
ucl3 <- rbind.data.frame(cbind.data.frame(setting=plan_names[1],result_data_geem_basic_ambulatory),
                 cbind.data.frame(setting=plan_names[2],result_data_geem_GP_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[3],result_data_geem_GP_residential),
                 cbind.data.frame(setting=plan_names[4],result_data_geem_WO_intensive_ambulatory),
                 cbind.data.frame(setting=plan_names[5],result_data_geem_WO_residential)
) %>% 
    dplyr::group_by(setting) %>% 
    dplyr::filter(model=="No time-varying coefficients, lag=1") %>% pull(`CI.upper`) %>% as.numeric()


# Calcula los logaritmos de los RRR y de los límites de sus intervalos de confianza
log_rrr3 <- log(rrr3)
log_lcl3 <- log(lcl3)
log_ucl3 <- log(ucl3)

# Calcula las varianzas
variances3 <- ((log_ucl3 - log_rrr3) / qnorm(0.975))^2

# Meta-análisis usando un modelo de efectos aleatorios
meta_analysis3 <- rma(yi=log_rrr3, sei=sqrt(variances3), method="REML")

# Prueba de Cochran Q para evaluar la heterogeneidad
cochrans_q3 <- meta_analysis3$QE
p_value_q3 <- meta_analysis3$QEp

cat(sprintf("No time-varying coefficients, lag=1- Cochran's Q = %.2f, p-value = %.4f\n", cochrans_q3, p_value_q3))
No time-varying coefficients, lag=1- Cochran's Q = 14.24, p-value = 0.0066
Code
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
###_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
###

# Comparaciones de Altman para todas las combinaciones de RRR
results <- expand.grid(i=1:5, j=1:5) %>%
  filter(i < j) %>%
  mutate(
    drrr = log_rrr[i] - log_rrr[j],
    var_drrr = variances[i] + variances[j],
    se_drrr = sqrt(var_drrr),
    z_value = drrr / se_drrr,
    p_value = 2 * (1 - pnorm(abs(z_value))),
    ratio_estimates = exp(drrr),
    lower_bound = exp(drrr - qnorm(0.975) * se_drrr),
    upper_bound = exp(drrr + qnorm(0.975) * se_drrr)
  ) %>% 
  dplyr::mutate(i= dplyr::case_when(i==1~plan_names[1], i==2~plan_names[2], i==3~plan_names[3], i==4~plan_names[4], i==5~plan_names[5]), j= dplyr::case_when(j==1~plan_names[1], j==2~plan_names[2], j==3~plan_names[3], j==4~plan_names[4], j==5~plan_names[5]))

results2 <- expand.grid(i=1:5, j=1:5) %>%
  filter(i < j) %>%
  mutate(
    drrr = log_rrr2[i] - log_rrr2[j],
    var_drrr = variances2[i] + variances2[j],
    se_drrr = sqrt(var_drrr),
    z_value = drrr / se_drrr,
    p_value = 2 * (1 - pnorm(abs(z_value))),
    ratio_estimates = exp(drrr),
    lower_bound = exp(drrr - qnorm(0.975) * se_drrr),
    upper_bound = exp(drrr + qnorm(0.975) * se_drrr)
  ) %>% 
  dplyr::mutate(i= dplyr::case_when(i==1~plan_names[1], i==2~plan_names[2], i==3~plan_names[3], i==4~plan_names[4], i==5~plan_names[5]), j= dplyr::case_when(j==1~plan_names[1], j==2~plan_names[2], j==3~plan_names[3], j==4~plan_names[4], j==5~plan_names[5]))

results3 <- expand.grid(i=1:5, j=1:5) %>%
  filter(i < j) %>%
  mutate(
    drrr = log_rrr3[i] - log_rrr3[j],
    var_drrr = variances3[i] + variances3[j],
    se_drrr = sqrt(var_drrr),
    z_value = drrr / se_drrr,
    p_value = 2 * (1 - pnorm(abs(z_value))),
    ratio_estimates = exp(drrr),
    lower_bound = exp(drrr - qnorm(0.975) * se_drrr),
    upper_bound = exp(drrr + qnorm(0.975) * se_drrr)
  ) %>% 
  dplyr::mutate(i= dplyr::case_when(i==1~plan_names[1], i==2~plan_names[2], i==3~plan_names[3], i==4~plan_names[4], i==5~plan_names[5]), j= dplyr::case_when(j==1~plan_names[1], j==2~plan_names[2], j==3~plan_names[3], j==4~plan_names[4], j==5~plan_names[5]))

# Imprime los resultados
#el 5 tiene el mayor efecto en omparación con el resto
#
bind_rows(cbind.data.frame(model= "No time-varying coefficients, no weight",results), 
                 cbind.data.frame(model= "No time-varying coefficients, lag=0",results2), 
                 cbind.data.frame(model= "No time-varying coefficients, lag=1",results3)
                 ) %>% 
  dplyr::mutate(across(c("drrr","var_drrr", "se_drrr", "z_value", "ratio_estimates", "lower_bound", "upper_bound"),~sprintf("%1.2f",.))) %>% 
  dplyr::mutate(across(c("p_value"),~sprintf("%1.4f",.))) %>% 
  dplyr::mutate(`rrr (95% IC)`= paste0(ratio_estimates, " (", lower_bound, ", ", upper_bound,")")) %>% 
  dplyr::select(-any_of(c("drrr", "var_drrr", "se_drrr", "z_value", "ratio_estimates", "lower_bound", "upper_bound"))) %>% 
  {
  #copiar_nombres2()
    write.table(., file = paste0(folder_path,"rr_contrasts_240625.csv"), dec=",", sep="\t")
    knitr::kable(., size=10, format="html", caption="GEE Models, Heterogeneity & contrasts by setting") %>% 
     kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>% 
     kableExtra::scroll_box(width = "100%", height = "375px")  
  }
GEE Models, Heterogeneity & contrasts by setting
model i j p_value rrr (95% IC)
No time-varying coefficients, no weight basic ambulatory GP intensive ambulatory 0.5812 0.99 (0.96, 1.02)
No time-varying coefficients, no weight basic ambulatory GP residential 0.0288 1.06 (1.01, 1.12)
No time-varying coefficients, no weight GP intensive ambulatory GP residential 0.0180 1.07 (1.01, 1.14)
No time-varying coefficients, no weight basic ambulatory WO intensive ambulatory 0.2098 1.04 (0.98, 1.11)
No time-varying coefficients, no weight GP intensive ambulatory WO intensive ambulatory 0.1395 1.05 (0.98, 1.12)
No time-varying coefficients, no weight GP residential WO intensive ambulatory 0.6052 0.98 (0.91, 1.06)
No time-varying coefficients, no weight basic ambulatory WO residential 0.0112 0.90 (0.84, 0.98)
No time-varying coefficients, no weight GP intensive ambulatory WO residential 0.0266 0.91 (0.84, 0.99)
No time-varying coefficients, no weight GP residential WO residential 0.0005 0.85 (0.78, 0.93)
No time-varying coefficients, no weight WO intensive ambulatory WO residential 0.0040 0.87 (0.79, 0.96)
No time-varying coefficients, lag=0 basic ambulatory GP intensive ambulatory 0.4239 0.98 (0.94, 1.03)
No time-varying coefficients, lag=0 basic ambulatory GP residential 0.0895 1.05 (0.99, 1.11)
No time-varying coefficients, lag=0 GP intensive ambulatory GP residential 0.0298 1.07 (1.01, 1.14)
No time-varying coefficients, lag=0 basic ambulatory WO intensive ambulatory 0.4805 1.03 (0.95, 1.12)
No time-varying coefficients, lag=0 GP intensive ambulatory WO intensive ambulatory 0.2636 1.05 (0.96, 1.15)
No time-varying coefficients, lag=0 GP residential WO intensive ambulatory 0.6656 0.98 (0.89, 1.07)
No time-varying coefficients, lag=0 basic ambulatory WO residential 0.0077 0.89 (0.81, 0.97)
No time-varying coefficients, lag=0 GP intensive ambulatory WO residential 0.0313 0.90 (0.83, 0.99)
No time-varying coefficients, lag=0 GP residential WO residential 0.0006 0.84 (0.77, 0.93)
No time-varying coefficients, lag=0 WO intensive ambulatory WO residential 0.0100 0.86 (0.77, 0.96)
No time-varying coefficients, lag=1 basic ambulatory GP intensive ambulatory 0.3486 0.98 (0.94, 1.02)
No time-varying coefficients, lag=1 basic ambulatory GP residential 0.0397 1.07 (1.00, 1.15)
No time-varying coefficients, lag=1 GP intensive ambulatory GP residential 0.0086 1.09 (1.02, 1.17)
No time-varying coefficients, lag=1 basic ambulatory WO intensive ambulatory 0.4305 1.03 (0.96, 1.11)
No time-varying coefficients, lag=1 GP intensive ambulatory WO intensive ambulatory 0.1919 1.05 (0.98, 1.13)
No time-varying coefficients, lag=1 GP residential WO intensive ambulatory 0.3783 0.96 (0.88, 1.05)
No time-varying coefficients, lag=1 basic ambulatory WO residential 0.0143 0.90 (0.83, 0.98)
No time-varying coefficients, lag=1 GP intensive ambulatory WO residential 0.0466 0.92 (0.85, 1.00)
No time-varying coefficients, lag=1 GP residential WO residential 0.0005 0.84 (0.76, 0.93)
No time-varying coefficients, lag=1 WO intensive ambulatory WO residential 0.0116 0.88 (0.79, 0.97)

Time for this code chunk to run: 0 minutes

Session info

Code
message(paste0("R library: ", Sys.getenv("R_LIBS_USER")))

R library: E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/renv/library/R-4.1/x86_64-w64-mingw32

Code
message(paste0("Date: ",withr::with_locale(new = c('LC_TIME' = 'C'), code =Sys.time())))

Date: 2025-02-16 23:53:28

Code
message(paste0("Editor context: ", path))

Editor context: C:/Users/CISS Fondecyt/Mi unidad/Alvacast/SISTRAT 2022 (github)/env

Time for this code chunk to run: 0 minutes

Code
sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
) %>% 
 knitr::kable(caption = "R packages", format = "html",
      col.names = c("Row number", "Package", "Version"),
    row.names = FALSE,
      align = c("c", "l", "r")) %>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>% 
  kableExtra::scroll_box(width = "100%", height = "375px")  
R packages
Row number Package Version
assertthat 0.2.1 CRAN (R 4.1.2)
backports 1.4.1 CRAN (R 4.1.2)
base64enc 0.1-3 CRAN (R 4.1.1)
BiocManager 1.30.18 CRAN (R 4.1.3)
biostat3 0.1.8 CRAN (R 4.1.3)
boot 1.3-28.1 CRAN (R 4.1.3)
broom 1.0.1 CRAN (R 4.1.3)
cachem 1.0.6 CRAN (R 4.1.2)
callr 3.7.2 CRAN (R 4.1.3)
cellranger 1.1.0 CRAN (R 4.1.2)
checkmate 2.1.0 CRAN (R 4.1.3)
chisq.posthoc.test 0.1.3 Github (ebbertd/chisq.posthoc.test@186d2ca6bbdba9fc19601aff4696ae1b85e7e0b0)
class 7.3-20 CRAN (R 4.1.3)
cli 3.4.1 CRAN (R 4.1.3)
clipr 0.8.0 CRAN (R 4.1.2)
clisymbols 1.2.0 CRAN (R 4.1.3)
cluster 2.1.4 CRAN (R 4.1.3)
coda 0.19-4 CRAN (R 4.1.2)
codetools 0.2-19 CRAN (R 4.1.3)
colorspace 2.0-3 CRAN (R 4.1.2)
CompQuadForm 1.4.3 CRAN (R 4.1.1)
cowplot 1.1.1 CRAN (R 4.1.2)
crayon 1.5.2 CRAN (R 4.1.3)
data.table 1.14.2 CRAN (R 4.1.2)
data.tree 1.0.0 CRAN (R 4.1.3)
DBI 1.1.3 CRAN (R 4.1.3)
dbplyr 2.2.1 CRAN (R 4.1.3)
deldir 1.0-6 CRAN (R 4.1.1)
DescTools 0.99.48 CRAN (R 4.1.3)
devtools 2.4.4 CRAN (R 4.1.3)
DiagrammeR 1.0.9 CRAN (R 4.1.3)
digest 0.6.29 CRAN (R 4.1.3)
dplyr 1.0.10 CRAN (R 4.1.3)
e1071 1.7-11 CRAN (R 4.1.3)
easyalluvial 0.3.2 CRAN (R 4.1.2)
ellipsis 0.3.2 CRAN (R 4.1.2)
emmeans 1.8.1-1 CRAN (R 4.1.3)
estimability 1.4.1 CRAN (R 4.1.3)
evaluate 0.16 CRAN (R 4.1.3)
Exact 3.2 CRAN (R 4.1.3)
expm 1.0-0 CRAN (R 4.1.2)
extrafont 0.18 CRAN (R 4.1.3)
extrafontdb 1.0 CRAN (R 4.1.1)
fansi 1.0.3 CRAN (R 4.1.3)
farver 2.1.1 CRAN (R 4.1.3)
fastmap 1.1.0 CRAN (R 4.1.3)
forcats 0.5.2 CRAN (R 4.1.3)
foreign 0.8-83 CRAN (R 4.1.3)
Formula 1.2-4 CRAN (R 4.1.1)
fs 1.5.2 CRAN (R 4.1.3)
future 1.28.0 CRAN (R 4.1.3)
future.apply 1.9.1 CRAN (R 4.1.3)
gargle 1.2.1 CRAN (R 4.1.3)
geeasy 0.1.2 CRAN (R 4.1.2)
geeM 0.10.1 CRAN (R 4.1.3)
geepack 1.3.9 CRAN (R 4.1.3)
generics 0.1.3 CRAN (R 4.1.3)
ggalluvial 0.12.5 CRAN (R 4.1.3)
ggforce 0.4.1 CRAN (R 4.1.3)
ggformula 0.10.2 CRAN (R 4.1.3)
ggplot2 3.4.1 CRAN (R 4.1.3)
ggridges 0.5.4 CRAN (R 4.1.3)
ggstance 0.3.6 CRAN (R 4.1.3)
gld 2.6.5 CRAN (R 4.1.3)
globals 0.16.1 CRAN (R 4.1.3)
glue 1.6.2 CRAN (R 4.1.3)
googledrive 2.0.0 CRAN (R 4.1.2)
googlesheets4 1.0.1 CRAN (R 4.1.3)
gower 1.0.0 CRAN (R 4.1.2)
gridExtra 2.3 CRAN (R 4.1.2)
gtable 0.3.1 CRAN (R 4.1.3)
hardhat 1.2.0 CRAN (R 4.1.3)
haven 2.5.1 CRAN (R 4.1.3)
highr 0.9 CRAN (R 4.1.3)
Hmisc 4.7-1 CRAN (R 4.1.3)
hms 1.1.2 CRAN (R 4.1.3)
htmlTable 2.4.1 CRAN (R 4.1.3)
htmltools 0.5.3 CRAN (R 4.1.3)
htmlwidgets 1.5.4 CRAN (R 4.1.2)
httpuv 1.6.6 CRAN (R 4.1.3)
httr 1.4.4 CRAN (R 4.1.3)
interp 1.1-3 CRAN (R 4.1.3)
ipred 0.9-13 CRAN (R 4.1.3)
jpeg 0.1-9 CRAN (R 4.1.1)
jsonlite 1.8.2 CRAN (R 4.1.3)
kableExtra 1.3.4 CRAN (R 4.1.3)
knitr 1.40 CRAN (R 4.1.3)
labeling 0.4.2 CRAN (R 4.1.1)
labelled 2.10.0 CRAN (R 4.1.3)
later 1.3.0 CRAN (R 4.1.2)
lattice 0.20-45 CRAN (R 4.1.1)
latticeExtra 0.6-30 CRAN (R 4.1.3)
lava 1.6.10 CRAN (R 4.1.2)
lifecycle 1.0.3 CRAN (R 4.1.3)
listenv 0.8.0 CRAN (R 4.1.2)
lme4 1.1-30 CRAN (R 4.1.3)
lmom 3.2 CRAN (R 4.1.2)
lmtest 0.9-40 CRAN (R 4.1.3)
lubridate 1.8.0 CRAN (R 4.1.2)
magrittr 2.0.3 CRAN (R 4.1.3)
MASS 7.3-58.1 CRAN (R 4.1.3)
mathjaxr 1.6-0 CRAN (R 4.1.3)
Matrix 1.5-1 CRAN (R 4.1.3)
MatrixModels 0.5-1 CRAN (R 4.1.3)
memoise 2.0.1 CRAN (R 4.1.2)
MESS 0.5.12 CRAN (R 4.1.2)
meta 7.0-0 CRAN (R 4.1.2)
metadat 1.2-0 CRAN (R 4.1.2)
metafor 4.2-0 CRAN (R 4.1.2)
mice 3.14.0 CRAN (R 4.1.2)
mime 0.12 CRAN (R 4.1.1)
miniUI 0.1.1.1 CRAN (R 4.1.3)
minqa 1.2.4 CRAN (R 4.1.2)
mitools 2.4 CRAN (R 4.1.2)
modelr 0.1.9 CRAN (R 4.1.3)
mosaicCore 0.9.2.1 CRAN (R 4.1.3)
muhaz 1.2.6.4 CRAN (R 4.1.2)
multcomp 1.4-20 CRAN (R 4.1.3)
MuMIn 1.46.0 CRAN (R 4.1.2)
munsell 0.5.0 CRAN (R 4.1.2)
mvtnorm 1.1-3 CRAN (R 4.1.1)
nlme 3.1-159 CRAN (R 4.1.3)
nloptr 2.0.3 CRAN (R 4.1.3)
nnet 7.3-18 CRAN (R 4.1.3)
numDeriv 2016.8-1.1 CRAN (R 4.1.1)
pacman 0.5.1 CRAN (R 4.1.3)
parallelly 1.32.1 CRAN (R 4.1.3)
pillar 1.8.1 CRAN (R 4.1.3)
pkgbuild 1.3.1 CRAN (R 4.1.2)
pkgconfig 2.0.3 CRAN (R 4.1.2)
pkgload 1.3.0 CRAN (R 4.1.3)
png 0.1-7 CRAN (R 4.1.1)
polspline 1.1.20 CRAN (R 4.1.3)
polyclip 1.10-4 CRAN (R 4.1.3)
prettyunits 1.1.1 CRAN (R 4.1.2)
processx 3.7.0 CRAN (R 4.1.3)
prodlim 2019.11.13 CRAN (R 4.1.2)
profvis 0.3.7 CRAN (R 4.1.3)
progressr 0.11.0 CRAN (R 4.1.3)
promises 1.2.0.1 CRAN (R 4.1.2)
proxy 0.4-27 CRAN (R 4.1.3)
ps 1.7.1 CRAN (R 4.1.3)
purrr 0.3.5 CRAN (R 4.1.3)
quantreg 5.94 CRAN (R 4.1.3)
R6 2.5.1 CRAN (R 4.1.3)
ragg 1.2.3 CRAN (R 4.1.3)
RColorBrewer 1.1-3 CRAN (R 4.1.3)
Rcpp 1.0.9 CRAN (R 4.1.3)
readr 2.1.3 CRAN (R 4.1.3)
readxl 1.4.1 CRAN (R 4.1.3)
recipes 1.0.1 CRAN (R 4.1.3)
remotes 2.4.2 CRAN (R 4.1.3)
renv 1.0.1 CRAN (R 4.1.2)
reprex 2.0.2 CRAN (R 4.1.3)
rlang 1.0.6 CRAN (R 4.1.3)
rmarkdown 2.16 CRAN (R 4.1.3)
rms 6.3-0 CRAN (R 4.1.3)
rootSolve 1.8.2.4 CRAN (R 4.1.2)
rpart 4.1.16 CRAN (R 4.1.3)
rstudioapi 0.14 CRAN (R 4.1.3)
Rttf2pt1 1.3.11 CRAN (R 4.1.3)
rvest 1.0.3 CRAN (R 4.1.3)
sandwich 3.0-2 CRAN (R 4.1.3)
scales 1.2.1 CRAN (R 4.1.3)
sessioninfo 1.2.2 CRAN (R 4.1.2)
shiny 1.7.2 CRAN (R 4.1.3)
showtext 0.9-5 CRAN (R 4.1.3)
showtextdb 3.0 CRAN (R 4.1.3)
SparseM 1.81 CRAN (R 4.1.1)
stringi 1.7.6 CRAN (R 4.1.2)
stringr 1.4.1 CRAN (R 4.1.3)
survey 4.1-1 CRAN (R 4.1.2)
survival 3.4-0 CRAN (R 4.1.3)
svglite 2.1.0 CRAN (R 4.1.2)
sysfonts 0.8.8 CRAN (R 4.1.3)
systemfonts 1.0.4 CRAN (R 4.1.2)
tableone 0.13.2 CRAN (R 4.1.3)
textshaping 0.3.6 CRAN (R 4.1.3)
TH.data 1.1-1 CRAN (R 4.1.3)
tibble 3.1.8 CRAN (R 4.1.3)
tidylog 1.0.2 CRAN (R 4.1.3)
tidyr 1.2.1 CRAN (R 4.1.3)
tidyselect 1.2.0 CRAN (R 4.1.3)
tidyverse 1.3.2 CRAN (R 4.1.3)
timeDate 4021.106 CRAN (R 4.1.3)
tweenr 2.0.2 CRAN (R 4.1.3)
tzdb 0.3.0 CRAN (R 4.1.3)
urlchecker 1.0.1 CRAN (R 4.1.3)
usethis 2.1.6 CRAN (R 4.1.3)
utf8 1.2.2 CRAN (R 4.1.2)
vcd 1.4-10 CRAN (R 4.1.3)
vctrs 0.5.2 CRAN (R 4.1.3)
viridisLite 0.4.1 CRAN (R 4.1.3)
visNetwork 2.1.2 CRAN (R 4.1.3)
webshot 0.5.4 CRAN (R 4.1.3)
withr 2.5.0 CRAN (R 4.1.2)
xfun 0.33 CRAN (R 4.1.3)
xml2 1.3.3 CRAN (R 4.1.2)
xtable 1.8-4 CRAN (R 4.1.2)
yaml 2.3.6 CRAN (R 4.1.3)
zoo 1.8-11 CRAN (R 4.1.3)

Time for this code chunk to run: 0.4 minutes

Code
reticulate::py_list_packages()%>% 
 knitr::kable(caption = "Python packages", format = "html",
      col.names = c("Package", "Version", "Requirement"),
    row.names = FALSE,
      align = c("c", "l", "r", "r"))%>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>% 
  kableExtra::scroll_box(width = "100%", height = "375px")  
Python packages
Package Version Requirement
absl-py 2.1.0 absl-py==2.1.0
asttokens 2.4.1 asttokens==2.4.1
astunparse 1.6.3 astunparse==1.6.3
audioconverter 2.0.3 audioconverter==2.0.3
autograd 1.6.2 autograd==1.6.2
autograd-gamma 0.5.0 autograd-gamma==0.5.0
beautifulsoup4 4.12.3 beautifulsoup4==4.12.3
Brotli 1.1.0 Brotli==1.1.0
certifi 2023.11.17 certifi==2023.11.17
cffi 1.16.0 cffi==1.16.0
charset-normalizer 3.3.2 charset-normalizer==3.3.2
clarabel 0.9.0 clarabel==0.9.0
click 8.1.7 click==8.1.7
cloudpickle 3.0.0 cloudpickle==3.0.0
colorama 0.4.6 colorama==0.4.6
comm 0.2.1 comm==0.2.1
contourpy 1.2.0 contourpy==1.2.0
cvxopt 1.3.2 cvxopt==1.3.2
cvxpy 1.5.2 cvxpy==1.5.2
cycler 0.12.1 cycler==0.12.1
debugpy 1.8.0 debugpy==1.8.0
decorator 4.4.2 decorator==4.4.2
delete-chrome-history-py 0.1.8 delete-chrome-history-py==0.1.8
easyocr 1.7.1 easyocr==1.7.1
ecos 2.0.13 ecos==2.0.13
editdistance 0.8.1 editdistance==0.8.1
efficientnet 1.0.0 efficientnet==1.0.0
essential-generators 1.0 essential-generators==1.0
et-xmlfile 1.1.0 et-xmlfile==1.1.0
executing 2.0.1 executing==2.0.1
fancyimpute 0.7.0 fancyimpute==0.7.0
ffmpeg 1.4 ffmpeg==1.4
ffmpeg-python 0.2.0 ffmpeg-python==0.2.0
filedir 0.0.3 filedir==0.0.3
filelock 3.13.1 filelock==3.13.1
flatbuffers 24.3.25 flatbuffers==24.3.25
fonttools 4.47.2 fonttools==4.47.2
formulaic 1.0.1 formulaic==1.0.1
fsspec 2023.12.2 fsspec==2023.12.2
future 0.18.3 future==0.18.3
gast 0.6.0 gast==0.6.0
git-filter-repo 2.45.0 git-filter-repo==2.45.0
google-pasta 0.2.0 google-pasta==0.2.0
graphviz 0.20.3 graphviz==0.20.3
grpcio 1.65.4 grpcio==1.65.4
gTTS 2.5.1 gTTS==2.5.1
h5py 3.11.0 h5py==3.11.0
idna 3.6 idna==3.6
imageio 2.34.2 imageio==2.34.2
imageio-ffmpeg 0.5.1 imageio-ffmpeg==0.5.1
imgaug 0.4.0 imgaug==0.4.0
iniconfig 2.0.0 iniconfig==2.0.0
interface-meta 1.3.0 interface-meta==1.3.0
ipykernel 6.29.5 ipykernel==6.29.5
ipython 8.20.0 ipython==8.20.0
jedi 0.19.1 jedi==0.19.1
Jinja2 3.1.3 Jinja2==3.1.3
joblib 1.4.0 joblib==1.4.0
jupyter_client 8.6.0 jupyter_client==8.6.0
jupyter_core 5.7.1 jupyter_core==5.7.1
keras 3.4.1 keras==3.4.1
Keras-Applications 1.0.8 Keras-Applications==1.0.8
keras-ocr 0.9.3 keras-ocr==0.9.3
kiwisolver 1.4.5 kiwisolver==1.4.5
knnimpute 0.1.0 knnimpute==0.1.0
lazy_loader 0.4 lazy_loader==0.4
libclang 18.1.1 libclang==18.1.1
lifelines 0.28.0 lifelines==0.28.0
llvmlite 0.41.1 llvmlite==0.41.1
Markdown 3.6 Markdown==3.6
markdown-it-py 3.0.0 markdown-it-py==3.0.0
MarkupSafe 2.1.4 MarkupSafe==2.1.4
matplotlib 3.8.2 matplotlib==3.8.2
matplotlib-inline 0.1.6 matplotlib-inline==0.1.6
mdurl 0.1.2 mdurl==0.1.2
mido 1.3.3 mido==1.3.3
ml-dtypes 0.4.0 ml-dtypes==0.4.0
more-itertools 10.2.0 more-itertools==10.2.0
moviepy 1.0.3 moviepy==1.0.3
mpmath 1.3.0 mpmath==1.3.0
multipledispatch 1.0.0 multipledispatch==1.0.0
mutagen 1.47.0 mutagen==1.47.0
namex 0.0.8 namex==0.0.8
natsort 8.4.0 natsort==8.4.0
nest-asyncio 1.5.9 nest-asyncio==1.5.9
networkx 3.2.1 networkx==3.2.1
ninja 1.11.1.1 ninja==1.11.1.1
nose 1.3.7 nose==1.3.7
numba 0.58.1 numba==0.58.1
numexpr 2.10.0 numexpr==2.10.0
numpy 1.26.3 numpy==1.26.3
openai-whisper 20231117 openai-whisper==20231117
opencv-python 4.10.0.84 opencv-python==4.10.0.84
opencv-python-headless 4.10.0.84 opencv-python-headless==4.10.0.84
openpyxl 3.1.4 openpyxl==3.1.4
opt-einsum 3.3.0 opt-einsum==3.3.0
optree 0.12.1 optree==0.12.1
osqp 0.6.5 osqp==0.6.5
packaging 23.2 packaging==23.2
pandas 2.2.0 pandas==2.2.0
pandas-flavor 0.6.0 pandas-flavor==0.6.0
parso 0.8.3 parso==0.8.3
patsy 0.5.6 patsy==0.5.6
pillow 10.2.0 pillow==10.2.0
platformdirs 4.1.0 platformdirs==4.1.0
pluggy 1.5.0 pluggy==1.5.0
polars 1.9.0 polars==1.9.0
proglog 0.1.10 proglog==0.1.10
prompt-toolkit 3.0.43 prompt-toolkit==3.0.43
protobuf 4.25.4 protobuf==4.25.4
psutil 5.9.8 psutil==5.9.8
pure-eval 0.2.2 pure-eval==0.2.2
pyarrow 15.0.0 pyarrow==15.0.0
pyclipper 1.3.0.post5 pyclipper==1.3.0.post5
pycparser 2.22 pycparser==2.22
pycryptodomex 3.20.0 pycryptodomex==3.20.0
pydotplus 2.0.2 pydotplus==2.0.2
pydub 0.24.1 pydub==0.24.1
Pygments 2.17.2 Pygments==2.17.2
pyjanitor 0.26.0 pyjanitor==0.26.0
PyMuPDF 1.24.9 PyMuPDF==1.24.9
PyMuPDFb 1.24.9 PyMuPDFb==1.24.9
pyparsing 3.1.1 pyparsing==3.1.1
PyPDF2 3.0.1 PyPDF2==3.0.1
pyreadr 0.5.0 pyreadr==0.5.0
pytesseract 0.3.10 pytesseract==0.3.10
pytest 8.3.1 pytest==8.3.1
python-bidi 0.6.0 python-bidi==0.6.0
python-dateutil 2.8.2 python-dateutil==2.8.2
pytube 15.0.0 pytube==15.0.0
pytube3 9.6.4 pytube3==9.6.4
pytz 2023.3.post1 pytz==2023.3.post1
pywin32 306 pywin32==306
PyYAML 6.0.1 PyYAML==6.0.1
pyzmq 25.1.2 pyzmq==25.1.2
qdldl 0.1.7.post1 qdldl==0.1.7.post1
regex 2023.12.25 regex==2023.12.25
requests 2.32.3 requests==2.32.3
rich 13.7.1 rich==13.7.1
rpy2 3.5.16 rpy2==3.5.16
scikit-image 0.24.0 scikit-image==0.24.0
scikit-learn 1.3.2 scikit-learn==1.3.2
scikit-survival 0.22.2 scikit-survival==0.22.2
scipy 1.11.4 scipy==1.11.4
scs 3.2.6 scs==3.2.6
seaborn 0.13.2 seaborn==0.13.2
semantic-version 2.10.0 semantic-version==2.10.0
setuptools-rust 1.8.1 setuptools-rust==1.8.1
shapely 2.0.5 shapely==2.0.5
six 1.16.0 six==1.16.0
soupsieve 2.5 soupsieve==2.5
SpeechRecognition 3.10.1 SpeechRecognition==3.10.1
spyder-kernels 2.5.2 spyder-kernels==2.5.2
stack-data 0.6.3 stack-data==0.6.3
statsmodels 0.14.1 statsmodels==0.14.1
sympy 1.12 sympy==1.12
target 0.0.11 target==0.0.11
tensorboard 2.17.0 tensorboard==2.17.0
tensorboard-data-server 0.7.2 tensorboard-data-server==0.7.2
tensorflow 2.17.0 tensorflow==2.17.0
tensorflow-intel 2.17.0 tensorflow-intel==2.17.0
tensorflow-io-gcs-filesystem 0.31.0 tensorflow-io-gcs-filesystem==0.31.0
termcolor 2.4.0 termcolor==2.4.0
threadpoolctl 3.4.0 threadpoolctl==3.4.0
tifffile 2024.7.24 tifffile==2024.7.24
tiktoken 0.5.2 tiktoken==0.5.2
torch 2.4.0 torch==2.4.0
torchaudio 2.4.0 torchaudio==2.4.0
torchvision 0.19.0 torchvision==0.19.0
tornado 6.4 tornado==6.4
tqdm 4.66.1 tqdm==4.66.1
traitlets 5.14.1 traitlets==5.14.1
translator 0.0.9 translator==0.0.9
typing_extensions 4.9.0 typing_extensions==4.9.0
tzdata 2023.4 tzdata==2023.4
tzlocal 5.2 tzlocal==5.2
urllib3 2.1.0 urllib3==2.1.0
validators 0.33.0 validators==0.33.0
watchdog 3.0.0 watchdog==3.0.0
wcwidth 0.2.13 wcwidth==0.2.13
websockets 12.0 websockets==12.0
Werkzeug 3.0.3 Werkzeug==3.0.3
whisper 1.1.10 whisper==1.1.10
wrapt 1.16.0 wrapt==1.16.0
xarray 2024.1.1 xarray==2024.1.1
youtube-dl 2021.12.17 youtube-dl==2021.12.17
yt-dlp 2024.7.9 yt-dlp==2024.7.9

Time for this code chunk to run: 0 minutes

Output

Code
invisible("Adicionales en E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/para_pres_cum_mean.R")
folder_path <- ifelse(dir.exists("E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/"),
                      "E:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/",
                      "G:/Mi unidad/Alvacast/SISTRAT 2022 (github)/_proposal_grant/2023/")
save.image(paste0(folder_path,"an_grant_23_24_4.RData"))
invisible("ültima vs. estable")
#load(paste0(folder_path,"1_an_grant_23_24_4.RData"))

Time for this code chunk to run: 0.1 minutes